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ABSTRACT 

A detailed ID model of the stratocumulus-topped marine boundary layer is described. The model has three 
coupled components: a microphysics module that resolves the size distributions of aerosols and cloud droplets, 
a turbulence module that treats vertical mixing between layers, and a multiple wavelength radiative transfer 
module that calculates radiative heating rates and cloud optical properties 

The results of a 12-h model simulation reproduce reasonably well the bulk thermodynamics, microphysical 
properties, and radiative fluxes measured in an ~500-m thick, summertime marine stratocumulus cloud layer 
by Nicholls. However, in this case, the model predictions of turbulent duxes between the cloud and subcloud 
layers exceed the measurements. Results of model simulations are also compared to measurements of a marine 
stratus layer made under gale conditions and with measurements of a high, thin marine stratocumulus layer. The 
variations in cloud properties are generally reproduced by the model, although it underpredicts the entrainment 
of overlying air at cloud top under gale conditions. 

Sensitivities of the model results are explored. The vertical profile of cloud droplet concentration is sensitive 
to the lower size cutoff of the droplet size distribution due to the presence of unactivated haze particles in the 
lower region of the modeled cloud Increases in total droplet concentrations do not always produce less drizzle 
and more cloud water in the model The radius of the mean droplet volume does not correlate consistently with 
drizzle, but the effective droplet radius does. The greatest impacts on cloud properties predicted by the model 
are produced by halving the width of the size distribution of input condensation nuclei and by omitting the effect 
of cloud-top radiative cooling on the condensational growth of cloud droplets. The omission of infrared scattering 
produces noticeable changes in cloud properties. The collection efficiencies for droplets <30-pm radius, and 
the value of the accommodation coefficient for condensational droplet growth, have noticeable effects on cloud 
properties. The divergence of the horizontal wind also has a significant effect on a 12-h model simulation of 
cloud structure. 

Conclusions drawn from the model are tentative because of the limitations of the ID model framework. A 
principal simplification is that the model assumes horizontal homogeneity, and, therefore, does not resolve 
updrafts and downdrafts. Likely consequences of this simplification include overprediction of the growth of 
droplets by condensation in the upper region of the cloud, underprediction of droplet condensational growth in 
the lower region of the cloud, and underprediction of peak supersaturations. 


1. Introduction 

Low-lying marine stratiform clouds cover a third of 
the ocean surface and play an important role in the 
earth’s radiative heat balance (Warren et al. 1988). Be- 
cause cloud-top temperatures are similar to surface 
temperatures, the longwave radiative impact of these 
clouds on the global heat budget is minor; however, 
because they reflect much more sunlight than the un- 
derlying ocean surface, they strongly affect the global 
albedo. It has been estimated that the global cooling 
that would result from a 4% increase in areal coverage 
by marine stratocumulus clouds would offset the ex- 
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pected warming from a doubling of atmospheric carbon 
dioxide (Randall et al. 1984). Twomey et al. ( 1984) 
demonstrated that the reflectivity of water clouds with 
modest optical depths increases as the abundance of 
cloud condensation nuclei (CCN) increases. Albrecht 
(1989) argued that the fractional coverage of marine 
stratocumulus clouds is regulated by the drizzle rate, 
which regulates and is regulated by CCN abundance. 
Ackerman et al. ( 1993) showed through model simu- 
lations that marine stratocumulus cloud layers may 
limit their own lifetimes by depleting CCN; this is be- 
cause vertical mixing in the stratocumulus-topped ma- 
rine boundary layer depends on cloud optical depth, 
which is itself dependent on droplet concentrations and 
therefore CCN. Evidence that increases in CCN con- 
centrations can increase cloud albedo and decrease 
drizzle is provided by linear high-albedo cloud features 
seen in marine stratiform clouds; these features, which 
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are often hundreds of kilometers long, are known as 
“ship tracks’' because they are caused by ships (Con- 
over 1966; Coakley et al. 1987; Radke et al. 1989; King 
et al. 1993). These results and observations indicate 
that cloud microphysics plays an important role in de- 
termining the global albedo of marine stratiform 
clouds. 

Models of the stratocumulus-topped marine bound- 
ary layer range in complexity from a single mixed layer 
scheme (Lilly 1968) to 3D large-eddy simulations 
(Deardorff 1980; Moeng 1986). Between these ex- 
tremes are multilevel ensemble-averaged turbulence 
models, in which the combined effects of all eddy sizes 
are parameterized. Turbulence closure for models in 
the latter category include higher-ordered closure 
schemes in 1 D ( e.g., Chen and Cotton 1 987; Bougeault 
1985) and in 2D ( e.g., Moeng and Arakawa 1980). 
The E— e closure method of Duynkerke and Driedonks 
( 1987) predicts gradient transfer coefficients for tur- 
bulent fluxes through the turbulent kinetic energy ( E ) 
and its viscous dissipation rate (c)- Although Holtslag 
and Moeng ( 1991 ) have described how to parameterize 
countergradient diffusion in convective boundary lay- 
ers, the present model simply treats turbulent transport 
through the assumption of downgradient diffusion. 

In most models of marine stratocumulus, thermo- 
dynamics and cloud microphysics are reduced to the 
solution of conservation equations for entropy and total 
water, which are partitioned into their components by 
bulk condensation schemes. Chen and Cotton ( 1987) 
parameterized the production of drizzle ( which they 
found to play a significant role in the turbulent structure 
of the boundary layer). Nicholls ( 1987) investigated 
the production of drizzle through explicit cloud micro- 
physics modeling, but his turbulence model was pre- 
scribed rather than predictive, and the nucleation of 
cloud droplets was not treated. 

In this paper we present a model developed for in- 
vestigating the effects of cloud microphysics on the 
dynamics and structure of the stratocumulus-topped 
marine boundary layer. We have coupled an E- e tur- 
bulence mixing model with a size-resolved aerosol and 
cloud microphysics model (Toon et al. 1988) and a 
sophisticated radiative transfer code (Toon et al. 
1989a). Our study focuses on cloud microphysics and 
optical properties, at the expense of a highly simplified 
treatment of air motions. One motivation for develop- 
ing such a model is to see how well it can reproduce 
the real atmosphere, despite the inherent simplifica- 
tions. Another reason for developing a model with sim- 
plified treatment of air motions is to investigate aerosol- 
cloud interactions that evolve over timescales that are 
prohibitively long to investigate with more sophisti- 
cated models and current computational constraints. 

Following a description of the model its outputs are 
compared to observations from the North Sea described 
by Nicholls ( 1984) and Nicholls and Leighton ( 1986). 


We then discuss some of the assumptions and sensitiv- 
ities of the model. 

2. Model description 

To represent the stratocumulus-topped marine 
boundary layer, we have developed a ID, horizontally 
homogeneous, Eulerian model that has three coupled 
components: aerosol and cloud microphysics, turbulent 
mixing, and radiative transfer. These components and 
some numerical computational issues are described in 
this section. 


a . Aerosol and cloud microphysics 


The cloud microphysics model treats two types of 
particles: condensation nuclei (CN ) and water droplets. 
The CN are haze particles in stable equilibrium with 
the humidity in each layer of the model. The particle 
size distributions are represented by C(r, z, t), where 
Cdr is the mean number concentration of particles with 
radii between r and r + dr at height z and time t . The 
analytic form of the particle continuity equation that 
the model solves for each particle size and type is 


dC Cd(pw) d 
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In this equation, p(z> t) is the mean density of dry air, 
w(z) the prescribed mean vertical air velocity, u, (r, z) 
a particle fall velocity, K{z y t) a gradient-transfer co- 
efficient, o p a turbulent Prandtl number ( the ratio of the 
eddy diffusivities for momentum and particles), S n (r, 
z, t) a source term of particles, R n (r , z, /) a particle 
removal rate, £ r (r, z, t) a condensational growth rate, 
and K c (r , r\ z, t) a coalescence kernel (symbols are 
listed in appendix A). The second term on the left side 
of Eq. ( 1 ) represents the horizontal divergence that 
compensates for any change in air density due to ver- 
tical convergence. The third term on the left side of ( 1 ) 
represents the divergence of the vertical flux due to 
advection and sedimentation, and the fourth term is for 
the divergence of the vertical flux due to turbulent dif- 
fusion. 

On the right side of ( I ), the source term accounts 
for the creation of entirely new CN, while conversions 
between CN and droplets appear as source and removal 
terms. The third term on the right side of ( 1 ) represents 
the flux divergence in radius space due to condensa- 
tional growth and evaporation of water. The first inte- 
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gral represents the creation of particles of radius r due 
to collisions between smaller particles, and the second 
integral represents the loss of particles of radius r due 
to collisions with other particles. The collision integrals 
are summed over the particle classes as described by 
Toon etal. (1988). 

The particle size distributions for CN and droplets 
are each divided into 50 bins with geometrically in- 
creasing size, such that the particle volume doubles be- 
tween successive bins, resulting in a radius grid that 
spans from r mm = 0.005 pm to r inax = 500 pm. Within 
each droplet size bin, the model solves a continuity 
equation for the total volume concentration of dis- 
solved cloud condensation nuclei (CCN) to allow cal- 
culation of the equilibrium reduction in droplet vapor 
pressure (the solute effect). Keeping track of the vol- 
ume of dissolved CCN allows the model to conserve 
solute mass. Because particles are circulated in the ma- 
rine boundary layer between unactivated haze particles 
below cloud base and cloud droplets within the cloud, 
the evaporation of droplets to haze particles must be 
treated carefully. The total volume of dissolved CCN 
within a droplet size bin is not enough information to 
adequately treat the evaporation of droplets from a 
varying distribution of CCN sizes. Following the ideas 
of Turco et al. ( 1979), our model carries the second 
moment of the CCN volume distribution within each 
droplet size bin, from which the width of the dissolved 
CCN distribution is calculated. The coagulation ex- 
pressions for CCN volume and volume-squared differ 
from those for droplet number. This difference occurs 
because collisions between particles cause the number 
of droplets to decrease, while the total particle volume 
remains unchanged, and the total volume-squared in- 
creases. 

The model calculates the concentration of water va- 
por by treating vapor exchange with droplets, and ver- 
tical transport by turbulent diffusion and advection. De- 
fining G(z, t) as the mean mass concentration (units 
of g cm" 3 ) of water vapor at height z and time /, the 
analytic form of the vapor continuity equation that the 
model solves is 


OG 

dt 
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P dz 
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where p w is the density of liquid water, and a h the tur- 
bulent Prandtl number for heat ( assumed to be the same 
as that for vapor ; Stull 1 988 ) . The right side of Eq. ( 2 ) 
represents vapor exchange with droplets (water in haze 
particles is ignored for all but the radiative calcula- 
tions). Below we describe the analytic forms of the 
terms appearing in these equations and the boundary 
conditions. 


1 ) Condensational growth 

Droplet condensational growth is treated within each 
model layer using an average value of the supersatu- 
ration (calculated from the predicted values of water 
vapor concentration and temperature), in which the 
calculation averages over updrafts and downdrafts. Be- 
cause of this averaging process, the model does not 
represent the horizontal variability of real clouds, in 
which supersaturations are expected to be higher in up- 
drafts than in downdrafts. The growth equation that we 
use treats the effects on droplet temperature of radiative 
transfer. We use a form similar to that employed by 
Barkstrom ( 1978 ): 


dr _ gjfl vap (S + 1 ~ A k A s + gl^ratl) 

dt 1 + #«£lrt V ap 
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in which n vap is the saturation vapor pressure (ex- 
pressed as a number density, with units of mole- 
cules cm 3 ) , and S the supersaturation of water vapor 
(5 = \njn v9p \ - 1, where is the ambient vapor 
pressure: tu = GA v //V/ w ). The gas kinetic corrections 
to the diffusion coefficient of water vapor, D, and the 
thermal conductivity, K n are described by Toon et al. 
( 1989b). In these expressions the condensation coef- 
ficient (called the mass accommodation coefficient by 
Toon et al.) and the thermal accommodation coeffi- 
cients are taken to be unity. The factors F v and F y are 
ventilation corrections, which account for the effects of 
droplet sedimentation on the flux of water vapor mol 
ecules and heat to and from the droplets; we use the 
functions recommended by Pruppacher and Klett 
( 1978). The factors and A s in (3) account for, re- 
spectively, the increase in vapor pressure exerted by 
the droplet due to curvature (the Kelvin effect) and its 
reduction due to the solute effect. They are given by 




i ^ d ^ s m s M w 


(4) 


The solute composition determines the values of the 
remaining factors in A k ; is the dissociativity of the 
dissolved CCN (for ammonium bisulfate, = 2), <f> s 
the practical osmotic coefficient ( we use a value of 1 , 
which is valid for dilute solutions), and m s the mass of 
dissolved CCN, which is equal to the density of the 
solute times the average CCN volume in a droplet 
size bin. 

The q rMi term in Eq. (3) accounts for the radiative 
heating rate of droplets: 

</rad = J k abii (J — B)dlS y (5) 


where J(v) and B(v) are, respectively, the mean ra- 
diative intensity and the Planck function in a model 
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layer, £ ahs (^, r ) the absorption coefficient for a droplet 
(from Mie calculations, described below), and the in- 
tegration is over frequency, /a 

Treatment of the rapid time rate of change in the 
sizes of haze droplets would require exceedingly small 
time steps (small fractions of a second ), because their 
growth rates are so large in relation to their sizes. To 
allow the model to take longer time steps (typically 
on the order of seconds), the growth of haze droplets 
in subsaturated air is not treated explicitly. Thus, the 
Kelvin and solute factors are ignored when evaluating 
the growth rate for haze droplets (for which A k A s < 
I). By ignoring A k A s , all haze droplets evaporate 
when the relative humidity is <100%; ultimately they 
are transformed from droplets to CN, as described 
below. 

2) CCN ACTIVATION 

The CN in a size bin are activated to CCN when 5 
exceeds S cm , where the value of the critical supersat- 
uration ( S cn , ) corresponds to the mass of one CN par- 
ticle. In the absence of the radiative term </ rad , the max- 
imum value of S for which a droplet can maintain 
stable equilibrium with respect to condensation 
(S lt „,k ) can be calculated from a simplified form of the 
Kohler equilibrium relations ( Pruppacher and Klett 
1978). The corresponding droplet radius is r crjt . Be- 
cause we include the radiative term </ rad in the droplet 
condensation equation, we also take it into account in 
the determination of S crn . This is done by evaluating 
Sen, = ■S.-ri. K + £i£ 2 < 7 rad. where (he terms on the right 
side are evaluated at r CTll . A more rigorous treatment 
would be to determine a new r crU that also takes into 
account </ rad , but for the conditions and droplet sizes 
that we consider we have found this refinement to be 
unnecessary. Including q raii reduces the magnitude of 
S cn{ in regions of radiative cooling. Although 5 cri , can 
be <0, we do not allow CCN activation under sub- 
saturated conditions. 

Activation appears as a nearly instantaneous removal 
rate ( 10 3 s _1 ) in the CN equations, and it appears as a 
corresponding source term in the droplet and dissolved 
CCN equations. When a nucleated particle first appears 
in a droplet size bin it contains no water; however, the 
resulting large solute effect causes it to grow rapidly 
through its r cnt value. 

3 ) Total evaporation of droplets 

Condensation nuclei are created when the evapora- 
tion of a droplet to the next smallest size bin would 
remove all of the water. This occurs when the average 
volume of dissolved CCN in a droplet exceeds the 
droplet volume, or when the average volume-squared 
of dissolved CCN exceeds the square of the droplet 
volume. When evaporation from a droplet bin would 
produce either of these conditions, the number of drop- 


lets and the volume of CCN still undergo condensa- 
tional loss, but instead of being a source to the next 
smaller droplet size bin, they are a source of CN. The 
dissolved CCN is distributed over the CN size distri- 
bution in a similar manner to that described by Turco 
et al. ( 1979). The method fits the dissolved CCN vol- 
ume to a lognormal probability distribution function, 
and distributes the CCN volume over the CN grid in a 
manner that conserves both volume and number. Since 
it would be unrealistic to allow accumulation of CN 
that are too small to have ever been activated, particles 
are only evaporated into CN bins of sizes larger than 
the smallest CN that has been activated to a droplet 
over the entire model domain in the previous hour of 
simulation. 

4 ) Particle collisions 

Collisions include the thermal coagulation of parti- 
cles due to their Brownian motion and the gravitational 
collection of particles due to differences in their fall 
speeds. In Eq. ( 1 ), K c is the sum of the thermal coag- 
ulation kernel ( Ay) and the gravitational collection ker- 
nel ( K g ). Our treatment of thermal coagulation uses the 
kernel given by Fuchs ( 1964) for spheres. Turbulence 
and phoretic forces are ignored. For gravitational col- 
lection we use the standard definition for the geometric 
collection efficiency between particles of radius r and 
r': 

/C g (r, r') = 7r(r + r') 2 E(r , r\)|u f (r) - v f ( r')|, (6) 

where E(r, r ' ) is the geometric collection efficiency 
between the particles (equal to the product of the col- 
lision and coalescence efficiencies), and v f a particle 
fall speed. The collision efficiencies are interpolated 
from the values used by Hall (1980), which derive 
from a number of sources. For coalescence efficiencies 
we have used the formulation of Beard and Ochs 
(1984). 

5) Transport 

Particles are transported by three processes: sedi- 
mentation, turbulent diffusion, and advection. The 
expressions used for sedimentation velocities apply 
to two regimes. For Reynolds numbers <10 2 (r 
^ 10 gm) we employ a Stokes-Cunningham ex- 
pression (Toon et al. 1989b). For larger Reynolds 
numbers there are no analytic expressions for the fall 
speeds; we use the interpolation developed by Beard 
( 1976). For turbulent diffusion we set the turbulent 
Prandtl numbers for particles and heat to unity. Gra- 
dient transfer coefficients are taken from the turbu- 
lence model; their evaluation is described in the sec- 
tion on turbulent mixing. 

Advection is specified through a fixed profile of ver- 
tical wind ( w), which is calculated from the prescribed 
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divergence rate of the horizontal wind (div). Assuming 
that the divergence rate is independent of height, the 
vertical wind profile is 

div(/7 surf ~ p) 

w = ( / ) 

HPr-n 

[ Equation (7 ) simplifies to vv ~ div • z in the boundary 
layer]. Because we represent continuity equations in 
flux form, the air continuity equation appears implic- 
itly. Under conditions of vertical velocity convergence, 
air will accumulate if unbalanced by horizontal diver- 
gence, thereby altering the predicted state variables (in 
which p is implicit). To avoid an increase of air den- 
sity, the vertical flux convergence of air is balanced by 
a horizontal divergence term [the second term on the 
left side of Eq. ( 1 ) ] . In a sensitivity test presented be- 
low, air is allowed to accumulate (as it does when sur- 
face pressure builds ) by omitting this horizontal diver- 
gence term. 

6 ) New particle creation 

New sulfate particles are introduced to the model 
domain by specifying a production rate at each grid 
point, which is described by a lognormal size distri- 
bution of CN. Because the growth of embryonic CN 
due to gas-to-particle conversion is not treated, these 
particles are introduced at accumulation mode sizes (r 
— 0.1 ^m). The parameters that must be specified are 
the particle production rate, and the geometric mean 
number radius (r„) and standard deviation (a) of the 
CN distribution. The values chosen are described 
below. 

7 ) Boundary conditions 

The flux of particles across the lower boundary of 
the model is given by 

Fc = -C(Z|)[l/ d ep.p + to], (8) 

where Z\ is the altitude of the midpoint of the layer. 
The only deposition velocity (u dcpp ) that we consider 
in Eq. (8) is that due to gravity; phoretic forces are 
ignored. The deposition velocity is given by Giorgi 
( 1986), which treats transfer across the dynamic sub- 
layer where turbulent motions drive mixing, as well 
as transfer across the viscous layer where mixing is 
due to molecular motion. The deposition velocities de- 
pend upon the atmospheric stability and wind speed 
near the surface, and are linked to the surface bound- 
ary conditions of the turbulence model (described 
below). 

The flux of water vapor across the lower boundary 
of the model is given by 

F g = -Vdep.vlGUi) - GUs)]. (9) 

In this expression, v 6cp v is the deposition velocity for 
water vapor molecules, which we also evaluate using 


Giorgi’s ( 1986) parameterization. Giorgi’s parameter- 
ization involves several nondimensional numbers: for 
water vapor molecules both the drift velocity and the 
Stokes number are zero, the Schmidt number is the ra- 
tio between the diffusivilies of water vapor and air, and 
the heat transfer coefficient (described in appendix C) 
is used in place of the drag coefficient. Here G(Z|) is 
the concentration of water vapor in the lowest layer, 
and the concentration of water vapor at the surface, 
G(z s ), corresponds to 98.3% relative humidity at the 
temperature of the sea surface, in accordance with 
Raoulf s law. 

When the divergence rate is zero, particles and water 
vapor have no fluxes at the top of the model; otherwise 
the fluxes into the top of the model are set equal to the 
advective fluxes out of the bottom of the highest grid 
layer, so that the concentrations in the top layer are not 
changed by advection. 

b. Turbulent mixing 

1 ) Governing equations 

We use a gradient transfer approach to model tur- 
bulent fluxes. The evolution equations for the dynamics 
and thermodynamics of the cloud-topped boundary 
layer are expressed in flux form: 

d(pu) d(pw) d(wpu) d( du\ 

—-“~ar + ~sr - Jz r-'sj 

= pf(v-v f .) ( l()a) 

d(pv) d(pw ) d(wpv ) d ( dv\ 

= pf( w g - u) ( 10b) 

d(pO) 0(pw) O(pwO) d_ / Kn t?fl\ 

dt dz dz dz\o h ^ dz) 



The components of the geostrophic wind, and i/ g , 
are fixed, and independent of altitude. Thermodynam- 
ics is influenced by radiative transfer through the di- 
vergence of the net upward radiative flux. The effect 
of cloud microphysics appears through the phase 
change of G, the water vapor concentration. 

The essence of the E-t closure method is to diag- 
nose the values of K m from the predicted values of the 
turbulent kinetic energy E and its viscous dissipation 
rate t, through the relation K m = cJL 2 It. The values of 
c ^ and all the constants used in the turbulence model 
are taken from Duynkerke ( 1988). 

The turbulent kinetic energy equation used in the 
model is 
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(){pE) d(pw) djwpE) ^ 0 / K m 0E\ 

dt Oz Oz Oz \ o h ; ^ Oz ) 

' — , 

TRANS 


+ 



SHEAR 


+ w’0[ - pe. (11) 

“v0 


BUOY 


The sources of E are transport (TRANS), shear 
(SHEAR), and buoyancy (BUOY). The pressure-ve- 
locity covariance is implicit in the gradient transport 
component of TRANS ( Duynkerke and Driedonks 
1987). The buoyancy flux is computed from the tur- 
bulent flux of virtual potential temperature, which is 
defined by = 0( 1 + 0.6 1 q s , — q \ ) , where q v is specific 
humidity and q\ the specific liquid water content. For 
the temperate conditions we are considering, q x and q ] 
are indistinguishable from the vapor and liquid water 
mass mixing ratios, respectively: 


G _ G 
p + G p 


and 


Pi Pi 
q x = — — » - , 
Pi + P P 


( 12 ) 


where the concentration of liquid water in a layer is 
calculated by integrating the volume of water in the 
droplet size distribution. Because 0 V is not conserved 
under saturated conditions, gradient transport of i9 v is 
used only for unsaturated air. For saturated air it is 
assumed that perturbations in vapor correlate with 
perturbations in temperature according to the Clau- 
sius-Clapeyron relation, and an expression for 
BUOY can be derived from the gradient transfer of 
two semi-conservative variables: equivalent poten- 
tial temperature (0 t . ) and specific total water content 
(<y t ). Duynkerke and Driedonks ( 1987) and Duynk- 
erke ( 1988) used an “all or nothing' 1 condensation 
scheme for their buoyancy flux treatment; we use a 
partial condensation scheme as described by Bou- 
geault (1981 ). In the partial condensation scheme, 
the saturated and unsaturated fluxes are combined 
through a weighting factor that represents the frac- 
tion of saturated air at a given level. Our model as- 
sumes an exponential probability distribution of the 
total water and liquid potential temperature (0|), the 
covariance of which is taken from Mellor and Ya- 
mada ( 1982). Details for the evaluation of BUOY 
are given in appendix B. The partial condensation 
scheme can be used to diagnose the cloud liquid wa- 
ter in each model layer. However, because our model 
explicitly predicts cloud droplet size distributions, 
this diagnosed cloud liquid water is not used. 

The most difficult quantity to model is e. We use the 
flux form of the e equation given by Duynkerke 
(1988): 

d(pe) d(pw) , d(wpe) d ( K m de\ 

-^r~ i — + -i)r-Wz\v t p dz) 


= ~(ci <P - c 2t pe), (13) 
E 


in which c u and c 2i are constants, and P is the local 
production of E: 

P = SHEAR + max(0, BUOY) + max(0, TRANS). 

(14) 

The inclusion of the transport term in P was determined 
by Duynkerke ( 1988) to be a necessary extension to 
the standard E- e model, enabling the model to produce 
results from the neutral boundary layer in agreement 
with higher-order closure models. 

All of the governing equations are solved using the 
same numerical methods applied to the cloud micro- 
physics continuity equations described by Toon et al. 
( 1988), where the tracer concentration analogs in - co- 
ordinates are pw, pu, p0 , pE , and pe. Because E and e 
are evaluated at grid layer boundaries, their transport 
requires gradient transfer coefficients to be interpolated 
from their values at layer boundaries, for which we use 
linear interpolation. We also use linear interpolation 
between grid centers for other quantities needed for the 
evaluation of buoyancy fluxes at layer boundaries. Be- 
cause of the nonlinear dependence of saturation vapor 
pressure on temperature, we evaluate vapor pressure at 
the layer boundaries using the interpolated tempera- 
tures. 

2) Boundary conditions 

Surface momentum fluxes are evaluated through 
bulk transfer coefficients: 

E ru = -c d Mpu(zi) and F pv = - c ti Mpv(z\), ( 15) 

where c tl is a drag coefficient, «(zi) and u(<T) are the 
mean wind components, and M = [«(^i) 2 + t(^) 2 J ,/2 
is the mean wind speed. The evaluation of c d , which 
depends on the stability of the surface layer and the 
surface stress, is described in appendix C. For the sur- 
face heat flux we apply Giorgi’s ( 1986 ) formulation of 
deposition velocity to the transport of heat, analogous 
to our treatment of vapor: 

E p0 = - p ujcp.4 0( z .\ ) - 0( z, ) U (16 ) 

where z s refers to the surface value and the deposition 
velocity is nearly the same as that for vapor (except 
that the Schmidt number for heat is the ratio of the 
viscosity to the thermal diffusivity of air). 

The boundary values for E and e are calculated at an 
altitude of 1 .5 m (5% of the thickness of the lowermost 
grid layer) from the relations given by Duynkerke 
(1988): 



(17) 

(18) 
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in which the convective velocity scale vv* is zero under 
stable conditions, and under unstable conditions 

vv* = • (19) 

In the highest layer of the model, the values of w, v , 
and 0 are all fixed at their initial values. The fluxes of 
E and e are set to zero at the top of the model. 

c. Radiative transfer 

1 ) General description 

The radiative transfer model is used to calculate ra- 
diative heating rates for use in the thermodynamic 
equation and in the droplet condensational growth 
term, and to calculate cloud optical properties. Only an 
overview is presented here, since the details of the tech- 
niques are described by Toon et al. ( 1988). 

The model treats multiple scattering over 26 solar 
wavelengths (0.26 p m < X < 4.3 pun), and absorption 
and scattering over 14 infrared wavelengths (4.4 pm 
< X < 62 pm). Blackbody energy beyond those wave- 
length domains is included to agree with the Stefan - 
Boltzmann law (though neither absorption nor scatter- 
ing are treated beyond the wavelength domains). An 
exponential sum formulation is used to treat gaseous 
absorption coefficients. The optical properties of par- 
ticles are determined through Mie calculations, in 
which the complex refractive index for liquid water is 
used, as interpolated from the datasets of Painter et al. 
(1969), Palmer and Williams (1974), and Downing 
and Williams (1975). Mie calculations are averaged 
over six radius subdivisions for each particle size bin. 

The distributions of particles, water vapor, and tem- 
perature are taken from the cloud microphysics and tur- 
bulent transport modules. Carbon dioxide is specified 
to be 340 ppmv, and the ozone profile is taken from 
the U.S. Standard Atmosphere ( NOA A 1 976 ) . For the 
radiative transfer calculations, the CN distribution is 
transformed from dry radius to the wet radius in equi- 
librium with the average relative humidity in each 
layer, which is not allowed to exceed 99.9% for this 
calculation. 

2) Boundary conditions 

The grid for radiative transfer includes one layer that 
lies above the grid used in the other modules. This layer 
(indicated by 0) is specified to be devoid of particles, 
its temperature is taken to be the same as that in the 
underlying layer, and the column of water vapor in the 
extra layer is fixed. Because the isothermal profile is 
an unrealistic representation of the actual atmospheric 
structure, we impose a downwelling source of long- 
wave radiation into the top of layer 0. The solar zenith 
angle as a function of time is calculated using standard 
methods. The reflectivity of the surface is fixed at 7% 


for all solar wavelengths and incident angles; its emis- 
sivity is taken to be unity. 

d. Numerical issues 

For the sake of clarity we have presented the conti- 
nuity equations above using geometric altitude (z) as 
the vertical coordinate; operationally, however, the 
equations are cast and solved using nondimensional 
pressure (sigma) as the vertical coordinate. The nu- 
merical algorithms used by the model to solve the con- 
tinuity equations are described by Toon et al. ( 1988). 

Below the highest level of the radiative transfer 
model, the grid is divided into 50 layers. The layer 
nearest the surface is 30 m thick. Layer thickness de- 
creases upward to 10-m resolution at the initial altitude 
of the inversion, above which the thickness is constant 
(at 10 m). The increased resolution above the initial 
height of the inversion allows the boundary layer to 
deepen unimpeded by resistance that can be caused by 
decreasing the vertical resolution upward of the inver- 
sion. Using thinner layers did not affect the model re- 
sults. The variables on the grid are staggered as fol- 
lows: the concentrations of particles and vapor, as well 
as horizontal winds and temperature, are defined at 
layer midpoints; vertical wind, gradient transfer coef- 
ficients, turbulent kinetic energy and its dissipation 
rate, and radiative fluxes are defined at layer bounda- 
ries. 

To avoid numerical instabilities, a variable time step 
is used. At each time step, the concentration of droplets 
in bins within 1 % of the peak concentration on the grid 
is not allowed to change by more than a factor of two. 
Likewise, S is not allowed to change by more than a 
factor of two when its absolute value exceeds 10 4 . 
Finally, E is not allowed to change by more than a 
factor of two when greater than £ mjn = 10 10 m 2 s 2 
(£ is not allowed to decrease below £ min ). If variations 
exceed these limits, then all dynamic variables are reset 
to their values at the end of the previous time step, the 
length of the time step is reduced, and the calculation 
is repeated. In this manner the time step is allowed to 
vary between 5 X 10 * and 40 s. 

The exponential scheme used for advection on the 
spatial grid, described by Toon et al. (1988), is not 
used here because instabilities resulted from the advec- 
tion due to subsidence of sharp gradients at the inver- 
sion (where mixing is negligible). Instead, upwind ad- 
vection, which is numerically diffusive, is used. Be- 
cause the details of the entrainment of overlying air 
may be important to the long-term evolution of the stra- 
tocumul us-topped marine boundary layer, the diffusive 
advection scheme may affect model results over long 
simulations. 

3. Comparisons of model results with measurements 

To investigate the extent to which the model can 
reproduce observations, we used two sources of data. 
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The first of these is a detailed set of airborne measure- 
ments taken in a layer of stratocumulus cloud over the 
North Sea by Nicholls ( 1984) ( hereafter referred to as 
N84). These measurements include vertical profiles of 
cloud microstructure, bulk thermodynamics, radiative 
fluxes, and turbulence through the depth of the cloud- 
topped boundary layer. The second source of data, 
which is less detailed than that of N84, consists of two 
sets of airborne measurements taken in stratiform 
clouds over the North Sea by Nicholls and Leighton 
( 1986) (hereafter referred to as NL). 

a. Comparisons with the measurements of Nicholls 

(1984) 

l ) Method 

The airborne measurements described by N84 (also 
by NL as flight 526) were obtained on 22 July 1982. 
Because these measurements represent a snapshot in 
time (averaged between 1100 and 1450 LST) of a 
boundary layer air mass with an uncertain history, dif- 
ficulties are to be expected in reproducing them with a 
time-dependent model. In using this same dataset, 
Duynkerke and Driedonks (1987) initialized their 
model with the observed cloud structure and compared 
their model results with the measurements after 2 hours 
of simulation time, before the calculations departed too 
far from the initial state. Bougeault ( 1985) took an- 
other approach in comparing model results with similar 
observations: the model was initialized with a cloud- 
less, stable atmosphere, simulating 7 days of evolution 
(when steady state was reached), and the results on 
noon of the fourth day were compared with the mea- 
surements. We take an intermediate approach by ini- 
tializing our model with a cloudless, slightly unstable 
boundary layer at midnight and comparing the model 
results 12 hours later with the measurements of N84. 

Since N84 did not report on particles of radius < 1 
//m, we initialize the model with 1000 CN per cubic 
centimeter spread uniformly through the depth of the 
model, with sizes specified by a lognormal distribution 
with r n = 0.05 /im, and o = 2.5. The CN composition 
is assumed to be ammonium bisulfate (NH 4 HS0 4 ), in 
accord with Covert’s ( 1988 ) measurements of non -sea 
salt sulfate in the remote North Pacific. The CN pro- 
duction rate and initial concentration were chosen to 
attain steady droplet concentrations of the approximate 
number observed by N84. In Fig. 1 the CCN activation 
spectrum produced by this CN distribution is compared 
to those measured by Hudson and Frisbie ( 1991 ) under 
marine stratus over the Pacific; it can be seen that the 
slopes of the two spectra are generally similar. Match- 
ing the slope of an observed CCN activation spectra is 
equivalent to matching the shape of the observed CN 
size distribution (assuming the same solute composi- 
tion). At the lower supersaturations depicted in Fig. 1, 
the slope of the activation spectrum for our initial CN 
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Fig. 1 . Comparison of initial cumulative CCN activation spectrum 
used in the model calculations (solid line) with the measurements of 
Hudson and Frisbie (1991). The three sets of measurements were 
made on successive days, during which the total particle concentra- 
tion was decreasing with time. 


distribution is steeper than Hudson and Frisbie’ s mea- 
surements. The difference in slopes is consistent with 
adding a larger proportion of smaller particles when 
increasing the total particle concentration. The rela- 
tively greater increase of smaller particles at higher to- 
tal particle concentrations is to be expected because the 
dominant source (by number) of marine aerosol is the 
photochemical conversion of gaseous sulfur to small 
sulfate particles. To maintain particle concentrations 
over the 1 2-h simulation, CN are produced at a rate of 
0.018 cm" 1 s “ 1 in the boundary layer. 

Sea surface temperatures were not reported by N84, 
therefore, we fix it at the climatological July value for 
the measurement area of 288 K (Tucker and Barry 
1984). The temperature in the lowest layer of the 
model is initialized 1 K below the sea surface temper- 
ature, which induces a small surface buoyancy flux 
(initially 8 W m 2 ). The initial lapse rate follows the 
dry adiabat up to an inversion altitude of 800 m, where 
the temperature jumps to the observed 286 K, and is 
isothermal above. The relative humidity through the 
depth of the boundary layer is initially 98.3% (same as 
the surface value), and above the inversion the water 
vapor mixing ratio is a constant 5 g kg 1 ( the average 
of the observations). The wind profile is initially geo- 
strophic and independent of altitude, as reported by 
N84 (n g = 8.5 m s" 1 , =0). No subsidence rates are 

reported in N84; to limit upward entrainment of the 
boundary layer, we take div = 2.5 X 10 6 s" 1 for the 
simulation. This yields a subsidence rate of 0.2 cm s 1 
at the initial altitude of the inversion. We fix the surface 
pressure at 1013 mb (air does not accumulate in this 
model simulation ) . 

The model simulation begins at local midnight, the 
sun rises at 0400 local time and reaches a minimum 
zenith angle of 35° at noon. The temperature in layer 0 
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Fig. 2. Time evolution of <a) the optical depth of the modeled boundary layer (at wavelength ), (b) 

the liquid water mixing ratio (in g kg and (c) the turbulent kinetic energy (in units of 10 4 cm 2 s 2 ). 
Output is made every half-hour of simulation. 


( the layer in the radiative transfer scheme above the rest 
of the grid ) is fixed at 286 K, while the column of water 
vapor in layer 0 is fixed at the climatological value of 
2.9 g cm' 2 (midlatitude summer profile; Anderson et al. 
1986). Combined with the down welling blackbody flux 
at the top of the model, this structure yields a downward 
longwave flux of 270 W m 2 at 900-m altitude, in agree- 
ment with the radiative calculations of N84. While this 
is 15 W m“ 2 below the measured value, Nicholls and 
Leighton ( 1 986 ) state that there was a systematic error 
in that particular measurement, and that the value cal- 
culated by N84 is more reliable. 


2) Results 

The evolution of optical depth of the modeled 
boundary layer provides a bulk measure of the cloud 
layer that develops in the model simulation (Fig. 2a). 
Prior to cloud formation the optical depth is dominated 
by extinction from the haze particles in the moist 
boundary layer. After 1.5 h of simulation, enough va- 
por has diffused upward to saturate the air and con- 
dense into a cloud, as seen in the time -height contour 
plot of q x ( Fig. 2b ) . The release of latent heat “shocks’ ’ 
the modeled atmosphere, as a large release of buoyancy 
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Fig. 3. Comparison between measurements and modeled profiles of (a) T , (b) </ v , (c) and (d) 
The dotted line is the initialization and the solid line the model output at 12 h. The triangles 
are N84’s measurements at 1 100 LST and the diamonds horizontal flight-leg averages. 


drives the boundary layer beyond thermodynamic equi- sorbs solar radiation, longwave cooling from cloud top 
librium by transporting heat and vapor upward into the is offset, leading to a reduction in mixing between the 
cloud layer. The shock of cloud formation is seen in cloud and the subcloud layer, and a gradual decline in 
the turbulence, which peaks at 2 h (Fig. 2c); this is optical depth by 12 h (local noon), 
followed by a peak in optical depth —1 h later, and 

then a maximum in < 7 , ~ I h after that. The accumula- ( j ) Temperature, bulk water, and wind 
tion of <71 produces drizzle, which depletes the cloud 

layer of liquid water. As the drizzle evaporates below Figure 3 shows that during 12 h of evolution the 
cloud base, the subcloud layer stabilizes with respect thermodynamic profiles of the simulated cloud-topped 
to the cloud layer, thereby reducing the supply of vapor boundary layer evolve toward those measured by N84. 
and driving the boundary layer back below thermody- The differences between the model results and the mea- 
namic equilibrium. As the sun rises and the cloud ab- surements are generally within the measurement un- 
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certainties [the only reported estimates of uncertainties 
for thermodynamic measurements were made by Nich- 
olls and Leighton ( 1 986 ) for cloud-top jumps: ± 1 .0 K 
for potential temperature and ±0.5 g kg 1 for total wa- 
ter mixing ratio] . The simulated cloud layer has a peak 
q ] that is —20% less than the average cloud-top value 
reported for the N84 measurements by Nicholls and 
Leighton (0.53 g kg 1 ). The lapse rate throughout the 
boundary layer that is predicted by the model is less 
stable than the observations. Because supersaturations 
predicted by the model are only slightly below zero in 
the lower region of the cloud, the lower temperatures 
in the model produce less vapor above cloud base 
(measured to be 380 ± 80 m) than in the measure- 
ments, which enhances differences between modeled 
and measured 0 C values. Below cloud base, both q s and 
0 e are better mixed in the measurements than in the 
model results. The large variability above cloud top in 
the measured values of q v and is due to a strong 
vertical gradient in vapor. 

No wind profiles were given by N84, but the winds 
were reported as —8.5 ms -1 from the north at all lev- 
els. Wind variation across the inversion was reported 
as small and variable. In the model results, variations 
in horizontal wind components across the inversion are 
<1 m s ~ ! . The modeled winds are approximately 9.2 
ms 1 from the north through the cloud layer. Drag- 
induced shear near the surface results in a northerly 
component of 6.5 ms 1 and an easterly component of 
0.9 ms -1 at the midpoint of the lowest grid layer ( 15- 
m altitude). 

( ii ) Cloud microstructure 

Figures 4 and 5 show that the microphysical prop- 
erties of the simulated cloud layer at 12 h generally 
match the measurements. No estimates of statistical un- 
certainties for the microphysics measurements were re- 
ported by N84, although some estimates of instrumen- 
tal uncertainties and systematic errors were presented. 
The measured vertical distribution of the total number 
concentration of droplets ( N dr „ ps ) , taken from averages 
over flight profiles presented by Nicholls and Leighton 
(1986), is roughly constant with height, with an av- 
erage of 90 cm 3 (Fig. 4a). There is some ambiguity 
regarding the smallest measured particle sizes: in N84 
and Nicholls and Leighton ( 1986) the radius range is 
given as 1 -400 /xm, while in Nicholls ( 1987 ) the same 
particle data are referred to as having a lower radius 
cutoff of 2 fim. Our modeled profile corresponding to 
the total number of haze particles (CN transformed to 
their equilibrium sizes, as used in the radiative calcu- 
lations) and droplets larger than 0.9-/im radius matches 
the measurements, although there is a trend of decreas- 
ing number with altitude that is not present in the mea- 
surements. This trend with height is sensitive to the 
particle size cutoff, as evident in the dotted line in Fig. 
4a for which the radius cutoff is 2.2 //m: the number 


of droplets increases with height. For the larger particle 
size cutoff, the profile of droplet concentration corre- 
lates with the supersaturation profile in the cloud layer 
(Fig. 4c). While it is often observed that in marine 
stratiform clouds the number of droplets is constant 
with height (e.g., Slingo et al. 1982; Noonkester 1984; 
Nicholls and Leighton 1986), in our model results the 
trend is determined by the profile of large, unactivated 
haze particles (but our results are affected by a super- 
saturation profile that may be unrealistic, as discussed 
below ) . Noonkester ( 1 984 ) used a lower droplet radius 
cutoff of 0.23 /xm, thereby ensuring that unactivated 
haze particles were counted. The prediction by our 
model of large, unactivated haze particles above cloud 
base is supported by the measurements of Akagawa and 
Okada (1993), who found high concentrations of small 
droplets (r < 2.5 //m) near the base of stratus clouds. 
However, the height dependence of droplet concentra- 
tion and its sensitivity to the particle size cutoff may 
also be artifacts of our ID model. 

Average droplet sizes are reported by Nicholls and 
Leighton ( 1 986 ) through the radius of the mean droplet 
volume, r v = 1 3^/(4^^)] 1/3 . Because it is not 
weighted by the size of the droplets, r v is sensitive to 
the number of small droplets (unlike the volume- and 
area- weighted mean radii). Because was ob- 

served to be roughly constant, and q\ increases with 
height in cloud, r v was observed to increase with height 
(Fig. 4b). An increase of droplet size with height can 
be derived from updraft parcel models in which only 
growth by condensation is treated (e.g., chapter 13 of 
Pruppacher and Klett 1978). Our model, which rep- 
resents averages over both updrafts and downdrafts and 
treats droplet coalescence, closely matches the mea- 
sured profile of f v when the radius cutoff is taken to be 
0.9 fi m. The agreement is not as good when the radius 
cutoff is 2.2 ^m, although the tendency with height 
remains. 

The modeled profile of supersaturation (5) at 12 h 
(Fig. 4c) is peaked at cloud top, where it has a value 
of 0.037%. This peak supersaturation is lower than pre- 
vious estimates for stratiform clouds, which are —0.2% 
(Hudson 1983); this is consistent with the expected 
underestimate of peak supersaturation due to averaging 
over updrafts and downdrafts in our model ( see section 
3c). The increase of S with height in the cloud layer 
can be explained by considering the flux of (see 
below) from an Eulerian viewpoint. There are signifi- 
cant upward vapor fluxes due to turbulent diffusion at 
all levels within the boundary layer, but there is virtu- 
ally no turbulent transport across the inversion. There- 
fore, vapor tends to accumulate just below the inver- 
sion, where it produces a peak in S. Furthermore, the 
peak in radiative cooling of the air at cloud top (due to 
the maximum in the longwave radiative flux divergence 
at that altitude — see below ) also contributes to the 
peak S at cloud top. A consequence of this height de- 
pendence of S in the cloud layer is that smaller CN will 
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Fig. 4. Comparison between modeled profiles of (a) total number concentration of droplets 
(JVdmps) and (b) radius of mean droplet volume (r v ) at 12 h (dotted and solid lines are for radius 
cut offs of 2.2 and 0.9 //m, respectively) and profile-average measurements from Nicholls and 
Leighton ( 1986) (diamonds). The CN are included at their equilibrium sizes, as used in the radi- 
ative calculations, (c) Modeled supersaturation (S) profile at 12 h (note the scale change at S 
= 0). (d) Comparison between measured drizzle flux profile from N84 (diamonds) and model 
results (solid line is total, dash -dotted line is for droplets with r < 25 ^m, and dashed line is for 
r > 25 //m). 


nucleate droplets near cloud top rather than at cloud 
base, leading to increasing numbers of activated drop- 
lets with altitude (correlating with the profile of droplet 
concentration with a radius cutoff of 2.2 gm). Here S 
is slightly negative between 520 and 380 m altitude, 
below which it steadily diminishes to a value of nearly 
-0.07 in the lowest model layer, corresponding to a 


relative humidity of 93% (the relative humidity is fixed 
at 98.3% in the surface skin layer). 

The profile of average supersaturation within the 
cloud layer predicted by our model (which averages 
over updrafts and downdrafts) differs from that given 
by updraft parcel models, which typically predict peak 
supersaturations near cloud base. In Hall’s ( 1980) 2D 
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Fig. 5. Comparisons of measured droplet size distributions from horizontal flight averages in 
N84 (diamonds) and model results at 12 h (solid lines) at altitudes of (a) 730 m, (b) 480 m, (c) 
300 m. and (d) 90 m The CN are included at their equilibrium sizes, as used in the radiative 
calculations. 


model simulation (with explicit microphysics) of a 
warm, convective maritime cloud, maximum supersat- 
urations were predicted near cloud top, as well as near 
cloud base. Peak supersaturations near cloud top are 
also predicted by stratocumulus cloud modeling results 
from 3D large-eddy simulations with explicit micro- 
physics (Kogan et al. 1994). However, in contrast to 
our model results, Kogan’s stratocumulus model pre- 
dicts supersaturations ( in updrafts ) that peak near cloud 
base and decline upward (until the inversion is ap- 
proached). We will consider some consequences of the 
horizontal averaging done by our model in section 3c. 


Figure 5 compares the modeled size distributions of 
droplets and haze with horizontal flight leg averages 
measured by N84 at four different altitudes. Consid- 
ering the uncertainties in the measurements and the 
seven orders of magnitude in droplet concentrations, 
the model results and the measurements are in surpris- 
ingly good agreement at all altitudes. There are peaks 
in the modeled size distributions due to haze particles 
(r ~ 1 /mi) that did not fall within the measurement 
range. Near cloud top ( 730 m) the peak in droplet num- 
ber at r = 10 fim is not as pronounced in the model 
results as it is in the measurements. In the lower region 




15 Ami. 1995 


ACKKRMAN KT AL. 


1217 


of the cloud (480 m), the model underpredicts the num- 
ber of small droplets (r < 10 jjm). The model under- 
predicts the droplet concentrations in the tail of the dis- 
tribution (r > 100 pm) at all altitudes. At all altitudes 
above 90 m, the model overpredicts the number of 
droplets with radii between 20 and 60 //m . In another 
comparison with these measurements, the steady-state 
model of drizzle production used by Nicholls ( 1987) 
also consistently overpredicted the concentrations of 
droplets with radii between 25 and 60 pm. However, 
as noted by N84, droplets of radius between 20 and 60 
pm were severely undercounted in the measurements, 
so the discrepancy may be the result of a measurement 
error. 

The drizzle flux, defined by 

4tt 

drizzle = *J" Pw J v f (r')C(r')r r 'dr' (20) 

and calculated by N84 from the measured droplet size 
distributions shown in Fig. 5, is compared to that pre- 
dicted by our model in Fig. 4d. The modeled flux is up 
to 50% larger than that measured, but the measured 
fluxes are only accurate to ±40% (Turton and Nicholls 
1987). The flux falls off rapidly below the supersatu- 
rated region, resulting in fluxes smaller than the mea- 
surements below 300-m altitude. For the model results 
shown in Fig. 4d, the drizzle flux is broken down into 
contributions from droplets of radii smaller and larger 
than 25 pm. As in the measurements (for which the 
drizzle flux by droplet size is not shown here), the driz- 
zle flux is dominated by droplets <25 p/m in radius in 
the top —200 m of the cloud layer, while at lower al- 
titudes larger drops dominate the drizzle flux. 

( iii ) Radiative transfer 

It can be seen from Fig. 6 that the radiative profiles 
of the modeled boundary layer at 12 h are within the 
uncertainties of the measurements (N84 estimated the 
measurement error of solar fluxes to be ±20 W m“ 2 ). 
Nicholls and Leighton measured an average broadband 
solar albedo above cloud top of 0.62 (they corrected 
all their solar fluxes to correspond to a zenith angle of 
35° ). Our model predicts a value of 0.61 at 1-km alti- 
tude. While this agreement is surprisingly good, the 
model overpredicts the amount of down welling solar 
flux below — 600-m altitude (Fig. 6a). This difference 
in flux is consistent with the underprediction of small 
droplets (r < 10 yum) in this region. The model slightly 
overpredicts the downwelling solar flux above the 
cloud; this may be due to underestimated columns of 
ozone and/or water vapor in the overlying model at- 
mosphere (climatological values were used for both). 

The “measured” longwave fluxes shown in Fig. 6b 
are actually theoretical calculations based on measured 
microphysical and thermodynamic profiles made by 
N84; Nicholls and Leighton ( 1986) considered these 
to be more reliable than the direct measurements due 


to a systematic error in the radiometer calibration above 
cloud. Furthermore, the radiometer for measuring up- 
welling longwave fluxes failed. The difference between 
the modeled and “measured" upwelling longwave flux 
near the surface could be due to a lower sea surface 
temperature assumed by Nicholls and Leighton than we 
used. 

The modeled boundary layer absorbs 53 W m 2 of 
solar radiation, and loses 65 W m 2 of longwave ra- 
diation, yielding a net loss of 12 W m~ 2 . For the mea- 
sured profiles shown in Fig. 6a, 80 W m 2 of solar ra- 
diation are absorbed, and 79 W m 2 of longwave ra- 
diation are lost between 90- and 930-m altitude, 
yielding a net gain of 1 W m 2 . The difference between 
our model predictions and the measurements is prin- 
cipally attributable to the underpredicted solar extinc- 
tion in the middle of the modeled cloud layer. A result 
of this difference in the net radiative budget of the 
boundary layer is that temperatures (and, hence, water 
vapor mixing ratios) in the cloud layer are below the 
observed values (Fig. 3). 

The modeled profiles of radiative heating are also 
shown in Fig. 6. Net radiative cooling at cloud top is 
limited to a region —50 m deep in the model, while 
flux profiles calculated by N84 indicate that cooling is 
restricted to a region 30 m deep. Net cooling at cloud 
top drives mixing by destabilizing the cloud layer itself, 
and also by destabilizing the cloud layer with respect 
to the subcloud layer. In our model results, solar heat- 
ing primarily offsets the longwave cooling at cloud top. 
As mentioned above, less solar heating was predicted 
by the model than was measured, and the difference 
occurs primarily in the lower region of the cloud layer. 
Because solar heating in this region can stabilize the 
cloud with respect to the subcloud layer, and thereby 
diminish mixing between the layers, the comparative 
lack of solar heating in our model allows greater mixing 
than was measured between the cloud and subcloud 
layers ( see below ) . 

(iv) Turbulence 

The measured turbulence profiles (for which no 
measurement uncertainties were reported) are gener- 
ally reproduced by the model results at 12 h (Figs. 7 
and 8). Shown in Fig. 7a is the measured profile of £, 
taken from 2.6-km segments of N84's 1 100 LST de- 
scent, which finished in the least convective part of the 
observational area; E was not reported for the horizon- 
tal leg averages, so it is not clear how representative 
these measurements are. Compared to the data, our 
model overpredicts the peak value of E by —100% in 
cloud and by —50% near the surface. The observations 
show a decrease downward from 400 m ( the observed 
cloud base during that descent) to a minimum at 200 
m, while the model predicts a decrease downward from 
700 m to a minimum at 300 m. As shown in the E 
budget (Fig. 7b), turbulence is driven by the buoyancy 
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(a) (b) 




(c) 



(d) 



-10 -8 -6-4-2 0 _ 2 

Longwave heating rate (K h ) 


Fig. 6. Measured (a) solar and {b) longwave radiative fluxes from N84: diamonds and triangles 
are upwelling and downwelling values, respectively. Solid lines are upwelling model results, and 
dotted lines are downwelling model results (at 12 h). (c) and (d) The modeled solar and longwave 
heating profiles, respectively. 


flux in cloud, where it is offset by dissipation and trans- 
port. Both shear and transport produce turbulence just 
below the cloud, while near the surface, turbulence is 
driven by shear and buoyancy. A measured budget of 
E that resembles the modeled budget is presented by 
N84 (not shown here). However, the measurements 
indicate a region of buoyant consumption of E below 
cloud base that is not represented in the model results. 
Comparison with N84 is complicated because, accord- 
ing to Duynkerke and Driedonks ( 1987), the measured 
transport term does not include pressure-velocity cor- 
relations. 


The buoyancy flux is the turbulent flux of 0 V , 
which is a variable that is not conserved under sat- 
urated conditions due to phase changes of water. As 
discussed in section 2, the buoyancy flux is calcu- 
lated by using a fractional condensation weighting 
factor to combine the flux from gradient transfer of 
0 V with the flux derived through a saturation condi- 
tion. This weighting factor, which represents the 
fraction of air that is saturated, is shown in Fig. 7c. 
It is seen to hover near 100% in the top ~200 m of 
the cloud layer, fall off gradually to ~60% just above 
300-m altitude, and drop off more rapidly below 
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Fractional condensation (% of air saturated) w'0^ (K cm s ) 

Fig. 7. (a) Comparison of measured E profile from N84 at 1 100 LST (diamonds) and model 
results at 12 h (solid line), (b) Modeled E budget (solid line is SHEAR, dotted line is BUOY, 
dashed line is TRANS, and dash -dotted tine is e). (c) Fractional condensation weighting factor, 
expressed as a percentage, (d) Comparison of measured turbulent buoyancy fluxes from N84 
horizontal flight averages (diamonds) with modeled profile at 12 h (line). 


300 m, where S starts to decrease dramatically ponents of the buoyancy flux; these are compared to 
(Fig. 4c). the modeled components in Fig. 8. The modeled tur- 

Figure 7d shows the buoyancy fluxes measured bulent fluxes in Fig. 8 are calculated by combining 
during horizontal flight legs and the values predicted the unsaturated fluxes (from gradient transfer) and 
by our model at 12 h. The model predicts a general the saturated fluxes (from Eqs. ( B3 )-( B5) given in 
trend with altitude that is consistent with the mea- appendix B) with the saturation weighting factor 
surements, but predicts positive values over too deep used for the buoyancy flux (see appendix B ). Near 
of a region, consistent with the larger values of E cloud top the sensible heat flux is overpredicted by 
predicted by the model. From the horizontal flight the model; in the lower half of the cloud layer the 
segments, N84 presented measurements of the com- fluxes of water are also overpredicted. 
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-1.0-0.5 0.0 a5 1.0_|.5 2.0 
w'qj (g kg cm s ) 


Fie; 8. Comparisons of measured turbulent fluxes of (a) 0, (b) <y v , and (c) q^ from horizontal flight averages 
in N84 (diamonds) with modeled profiles at 12 h (lines), 


( v ) Summary 

The results presented above show that some of the 
features of the stratocumulus-topped marine boundary 
layer predicted by the model are in agreement with the 
measurements of N84, but that some features are not 
as well reproduced by the model. The bulk thermody- 
namics of the boundary layer and the profile of the 
drizzle flux are generally reproduced by the model, as 
are the radiative fluxes in the upper region of the cloud. 
The model slightly underpredicts the concentration of 
10-jim radius droplets near cloud top, and it underpre- 
dicts the number of droplets < 10-/im radius in the 
lower region of the cloud. This lack of small droplets 
is responsible for the comparative lack of solar extinc- 
tion in the lower region of the cloud. A lack of solar 
absorption is consistent with slightly lower tempera- 
tures in the cloud layer and greater mixing with the 
subcloud layer in the model results than in the mea- 
surements. Further reasons for discrepancies between 
the measurements and the model results are considered 
in section 3c. 

b. Comparisons with measurements of Nic hoi Is and 
Leighton ( 1 986 ) 

We will now compare results from simulations using 
the model described in this paper with two sets of air- 
borne measurements taken in stratiform clouds over the 
North Sea by Nicholls and Leighton ( 1986 ), hereafter 
referred to as NL. The first set of measurements from 
NL (flight 564 made on 15 December 1982) provides 
a good test of the versatility of the model because the 
turbulence was quite different from that in flight 526, 
which we have compared to model results in section 


3a. In flight 564, gale conditions produced a cloud that 
nearly filled the depth of the boundary layer (670 m 
deep cloud in an 85()-m boundary layer). In contrast, 
in flight 526 cloud-top radiative cooling produced the 
buoyancy that drove mixing in a boundary layer that 
was only half filled by the cloud layer. In the second 
dataset from NL (flight 528 made on 29 July 1982), a 
thin cloud layer ( depth 1 90 m) capped a deep boundary 
layer ( 1260 m), resulting in a greatly reduced drizzle 
flux. Since NL presented significantly less detailed data 
than N84, the comparisons with model results will be 
less comprehensive than presented above. Our main 
objective is to show that the model responds to signif- 
icant variations in boundary conditions in a manner 
consistent with observations. 

1 ) Flight 564 ( shear-driven stratus ) 

In flight 564 the measured average wind speed was 
27 ms -1 . The cloud top was smooth and featureless; 
NL classified it as stratus, compared to stratocumulus 
for the other cases they studied that exhibited cellular 
patterns associated with convective motions. The dif- 
ferent appearances of the clouds reflects the different 
processes that produced the turbulent mixing. In mod- 
eling this dataset, Duynkerke and Driedonks ( 1988) 
found that the large drizzle fluxes produced by the thick 
cloud in flight 564 had a significant effect on the ver- 
tical structure of the cloud layer. 

Shown in Table 1 are the boundary and initial con- 
ditions used for the model simulations of flight 564. 
The divergence rate was increased in this model sim- 
ulation to prevent upward growth of the boundary layer 
due to increased shear-driven turbulence. Because the 
observations of flight 564 were made in December 
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Tabll I. Summary of (a) boundary conditions and (h) initial conditions for model simulations. Flight 526 corresponds to the base case 
simulation; flights 564 and 528 correspond to the additional model simulations presented in section 3b. In flight 564, a distribution of CN 
(representing sea salt) is fixed in the lowest model layer, as described in the text. 


Parameter 

526 

Flight 

564 

528 


(a) Boundary conditions 



Geostrophic wind speed (m s 1 ) 

8.5 

27 

3.5 

Divergence rate of horizontal wind velocity (s ! ) 

2.5 x 10 * 

11x10 h 

2 x 10 

Sea surface temperature (°C) 

15 

12.5 

16 

Temperature above inversion (°C) 

13 

9.5 

13 

Vapor mixing ratio above inversion (g kg ') 

5 

6.2 

2 

Downwelling longwave flux above inversion (W 

m 2 ) 270 

245 

250 

Solar zenith angle at noon (°) 

35 

75 

37 

Production rate of CN (cm ? s ') 

0.018 

0* 

0.022 


(b) Initial conditions 



Height of inversion (m) 

800 

9(H) 

1 250 

Temperature of lowest model layer (°C) 

14 

13 

17 

Initial concentration of CN (cm *) 

1000 

10 

12(H) 


(flight 526 was in July), the sea surface and the air 
above the inversion are colder, and the minimum solar 
zenith angle is much greater than for flight 526. Be- 
cause aerosols in the marine boundary layer during pe- 
riods of high wind consist mainly of sea salt (e.g., 
O’ Dowd and Smith 1993), and because droplet con- 
centrations measured during flight 526 generally de- 
creased with height (Fig. 9c), for this simulation we 
did not prescribe a distributed source of sulfate particles 
throughout the boundary layer, as we did in all the other 
model simulations. Instead, we fixed a size distribution 
of sea salt particles in the lowest layer of the model, 
based on measurements for 17 ms -1 wind speeds 
(O'Dowd and Smith 1993). The fixed distribution in 
the lowest layer provides CN to the boundary layer by 
mixing upward. For the distribution in the lowest layer, 
we used a total CN concentration of 200 cm \ a geo- 
metric mean number radius of 0.2 p m, and a geometric 
standard deviation of 2.5. We also initialized a sulfate 
distribution with a concentration of 10 cm * through- 
out the model domain with the same mean size and 
distribution width as in the base case simulation. Be- 
cause the dynamic relaxation time for such strong 
winds is about one day, we started the simulation at 
noon and ran the model for 27 h (measurements in 
flight 564 were between 1 100 and 1500 LST). 

In the 27-h model simulation, a cloud layer 600 m 
deep developed in an 800-m boundary layer (cloud 
depth in the model is taken to be the difference between 
the altitudes at which visibility is <1 km). Figure 9a 
shows that the model overpredicts the liquid water in 
the top 200 m of the cloud layer, suggesting that the 
model underpredicts cloud-top entrainment in this sim- 
ulation. The broadband solar albedo (above cloud top) 
of the modeled cloud at noon was 70%, compared to 
the measured value of 68%. As in the measurements, 


the vertical mixing in the model was driven primarily 
by shear, though the model predicted a region of pos- 
itive buoyancy flux below cloud top that was inconsis- 
tent with the measurements (not shown ). Although the 
modeled temperature profile matches the measured 
profile through most of the depth of the boundary layer 
(not shown), near cloud top the measured profile was 
more stable than indicated by the model. The over 
prediction of the instability of the temperature profile 
near cloud top indicates that the entrainment of inver- 
sion air is underpredicted by the model. The modeled 
instability is attributable to radiative cooling and results 
in a positive buoyancy flux. The shape of the modeled 
wind profile matches the measurements, although the 
observed jumps at cloud top are underpredicted by the 
model ( Fig. 9b ) . The increased wind shear at cloud top 
in the observations resulted in more entrainment of in- 
version air than in the model, which evaporated more 
droplets in the upper region of the cloud as well as more 
completely offsetting radiative cooling ( which was fur- 
ther reduced due to decreased liquid water). This in- 
terpretation is supported by comparing modeled and 
measured droplet concentrations ( Fig. 9c ) : the average 
number of droplets is captured well by the model, but 
there is more variation with height in the measure- 
ments. The local maximum in droplet concentration for 
droplets of radius >0.9 ^xm at ~350-m altitude in the 
model results is attributable to swollen haze particles 
(because of the high relative humidity); this is evident 
because there is no maximum in the concentration of 
droplets of radius >2.2 /im at that altitude. Finally, 
comparison of modeled and measured drizzle fluxes 
(Fig. 9d) shows that the modeled predictions exceed 
the measurements within the cloud layer. However, the 
differences are generally within the 40% uncertainty of 
the measurements (Nicholls 1987). 
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(c) (d) 


0 20 40 60_ 80 100 0.0 0.5 1.0 1.5 2.0 2.5 3.0 

N drops ( cm ) Drizzle flux (mm ) 

Fig. 9. Comparison between measurements (symbols) from flight 564 (Nicholls and Leighton 
1986) and modeled profiles (lines) at 25 h of (a) q x \ (b) wind speeds (dotted and solid lines are 
modeled u and v components, respectively; triangles and diamonds are the measured u and v 
components, respectively); (c) total number concentration of droplets (W (lf<)ps ) (dotted and solid 
lines are model output for radius cut offs of 2.2 and 0.9 //m, respectively); and (d) drizzle flux 
(solid line is model output). The measurements in (a), (b), and (d) are horizontal averages, the 
measurements in (c) profile averages. 

In summary, the model simulation for flight 564 re- thin cloud produced little drizzle. In the measured pro- 
produces the general features of the measurements: a files of T and q v (Figs. 10a and 10b, respectively) and 

deep cloud layer driven by wind shear, resulting in large the winds ( not shown), it is seen that the air below 700 

drizzle fluxes. However, entrainment at cloud top was m was strongly decoupled from the overlying air ( “de- 
less pronounced in the model than in the measurements, coupled 11 refers to a lack of vertical mixing between 

layers). Below 700 m the temperature profile was sta- 

2 ) Flight 528 ( high, thin stratocumulus) ble and the winds stronger. This structure reflects hor- 

izontal variations in air properties (the strong wind 

The cloud layer on flight 528 was measured to be shear implies different source regions for the air above 
190 m deep, atop a boundary layer 1260 m deep. This and below 700 m). Although time variations in the pro- 
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tiles of air properties seems likely, no information on 
time dependence was provided by NL. The structure of 
the observed boundary layer implies that there was little 
transport of water vapor from the sea surface to the air 
above 700 m. To restrict mixing (in the context of a 
gradient transfer formulation of turbulent fluxes) be- 
tween the upper and lower regions of the boundary over 
the 12-h model simulation, water vapor was initialized 
with a constant mixing ratio of 6.5 g kg 1 throughout 
the boundary layer (at the top of the boundary layer, 
where this much vapor would have saturated the air, 
we limited the initial relative humidity to 99%). Other 
initial and boundary conditions are listed in Table 1 . 

After 12 h of simulation, the model produces a cloud 
layer 210 m deep that caps a boundary layer 1200 m 
deep. As seen in Fig. 10a, the model underpredicts the 
temperatures in the upper region of the boundary layer; 
in the same region, the model overpredicts the water 
vapor mixing ratio (Fig. 10b). However, this apparent 
inconsistency is within the errors of the measurements. 
Neither the stable temperature profile near the surface 
nor the enhanced water vapor near the surface are re- 
produced by the model simulation. As mentioned pre- 
viously, these features of the profiles probably reflect 
horizontal variations (which cannot be represented in 
a ID model ) and/or unmeasured time variations in air 
properties. The model results closely match the mea- 
sured liquid water profile (Fig. 10c). There were 140 
droplets cm 3 measured in the cloud layer, and the 
model produced a cloud with 150 droplets cm \ The 
average broadband solar albedo measured above cloud 
top was 43%; the modeled value was 40%. Finally, as 
shown in Fig. lOd, the model simulation resulted in a 
drizzle flux that was much less than in the other model 
simulations ( Figs. 4d and 9d ), yet the measured drizzle 
flux for flight 528 (albeit only a single data point) was 
even smaller. 

In summary, the model simulations match the mea- 
surements of the thin cloud studied in flight 528. More 
structure was measured in the lower region of the 
boundary layer than the model indicated, but this dif- 
ference could have been due to horizontal variations 
and/or time dependence of the air properties. 

c. Likely consequences of ID model assumptions 

While the model results appear to generally match 
the measurements within the ranges of measurement 
uncertainties, it is worth considering the effects of some 
of the fundamental assumptions made in the model’s 
formulation. Foremost of these assumptions is that in 
our ID model, grid layers represent averages over up- 
drafts and downdrafts, whereas in a real cloud, the mi- 
crophysical environment of a downdraft differs from 
that in an updraft. To the extent that updrafts and down- 
drafts maintain their identities throughout deep regions 
of stratocumulus cloud layers (as in cumulus cloud 
fields), a ID model cannot represent the microphysical 


environment realistically. However, the ID model ap- 
proach is justifiable to the extent that stratocumulus 
cloud layers are horizontally homogeneous and the mi- 
crophysical properties are similar in updrafts and 
downdrafts. Thus, our model results represent the con- 
sequences of averaging over updrafts and downdrafts. 

Although there are no measurements of supersatu- 
ration profiles, numerical cloud models that consider 
the microphysical environment of updrafts ( such as up- 
draft parcel models and 2D and 3D models) predict 
peak supersaturations near cloud base within updraft 
regions. In contrast, our model does not predict a ( hor- 
izontally averaged) peak supersaturation near cloud 
base. One consequence is that the model probably un- 
derpredicts the condensational growth of some of the 
droplets in the lower region of the cloud layer (those 
in supersaturated updrafts, which are not represented 
in the model), resulting in an underprediction of small 
droplet concentrations in this region. This interpreta- 
tion of the discrepancies between the model results and 
measurements is consistent with the lack of droplets ( / 
> 2.2 pm) near cloud base in the model results in com- 
parison to the measurements (Fig. 4a). Furthermore, 
overprediction of evaporation in the lower region of 
cloud probably affects mixing, by destabilizing the 
lower region of the cloud with respect to the subcloud 
layer, thereby inducing a buoyancy flux that allows 
more mixing to occur than was observed, which is con- 
sistent with the comparisons with measured turbulent 
fluxes shown in Figs. 7d and 8. 

Although a peak in average supersaturation near 
cloud top may be physically representative of strato- 
cumulus cloud layers (Kogan 1994), it is likely that 
there are unsaturated regions in downdrafts near cloud 
top. In this case, by averaging over updrafts and down- 
drafts, the model neglects any time spent by droplets 
near cloud top in unsaturated downdrafts (much as it 
neglects time spent in saturated updrafts near cloud 
base). Hence, this model simplification probably re- 
sults in droplets being exposed to saturated conditions 
for artificially long durations near cloud top (much as 
they are exposed to artificially long durations of unsat- 
urated conditions near cloud base). This interpretation 
of likely model shortcomings is consistent with an un- 
derprediction of the peak in the droplet size distribution 
at r ~ 10 pm near cloud top (Fig. 5a). Artificially 
enhanced condensation near cloud top could result in 
overpredicting the drizzle fluxes there, by growing 
droplets into size ranges where they are more effi- 
ciently captured gravitationally by larger drops, con- 
sistent with overpredicted drizzle rates near cloud top 
in all the simulations (Figs. 4d, 9d, and lOd). In com- 
bination with an underpredicted drizzle flux below 
cloud base, the vertical gradient of drizzle flux is over- 
predicted, consistent with the comparison against ob- 
servations (Fig. 4d). 

Another likely consequence of averaging over up- 
drafts and downdrafts is that peak supersaturation val- 
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(a) (b) 




(c) (d) 




Fig. 10. Comparison between measurements (triangles) from flight 528 (Nicholls and Leighton 
1986) and modeled profiles at 12 h (lines) of (a) T, (b) </ v , (c) and (d) drizzle flux. The dotted 
lines in (a) and (b) are model initializations. The triangles in (a) -(c) are measurements from a 
profile at 1130 LST, the diamond in (d) is an average cloud-top measurement. 


ues are underpredicted by the model, due to the aver- 
aging process. Therefore, the model likely underpre- 
dicts the fraction of CN that are nucleated. To 
compensate for this shortcoming, the total number of 
CN input to the model is probably artificially high. 

We will now investigate the sensitivity of model out- 
puts to some input parameters and assumptions. 

4. Sensitivity of the model to physical assumptions 

Some of the physical processes represented in our 
model, such as the radiative effect on the condensa- 


tional growth of droplets and the scattering of radiation 
in the infrared, are not generally included in cloud mi- 
crophysical models. Therefore, we start this section by 
describing the effects of these processes on model re- 
sults. This is followed by investigating the sensitivity 
of the model outputs to the assumed values of the con- 
densational coefficient and the gravitational collection 
efficiencies, the flux-profile relationships used in the 
evaluation of surface boundary conditions, the proba- 
bility distribution function used in the partial conden- 
sation scheme, and the treatment of horizontal velocity 
divergence. For each of the sensitivity tests, the model 
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simulation of section 3a ( N84 and flight 526 in NL) is 
repeated with either a physical process omitted or a 
parameter changed. 

a . Radiative effects on condensational growth of 
droplets 

Several researchers have evaluated theoretically the 
effects of radiative heating on the condensational 
growth of droplets (e.g., Barkstrom 1978; Davies 
1985). Caughey and Kitchen (1984) used a parcel 
model to analyze measured droplet size distributions in 
a stratocumulus cloud layer and found that the effects 
of radiative transfer on droplet condensational growth 
were significant. However, these effects have been gen- 
erally ignored even in rather sophisticated models of 
cloud microphysics (e.g.. Hall 1980; Flossman et at. 
1985; Nicholls 1987). Near the top of a cloud layer, 
where radiative cooling rates are significant, cloud 
droplets are cooled relative to the ambient air by radi- 
ation. There are two direct effects of such temperature 
differences on cloud microphysics: increased conden- 
sational growth rates and a decrease in S cn( . These two 
processes act in opposite directions on droplet nuclea- 
tion: the first decreases the supersaturation, thereby de- 
creasing the number of droplets nucleated, while the 
second increases the number of droplets nucleated. 

To investigate the importance of this cloud droplet 
radiation effect, we ran the model without including 
the effects of radiative droplet cooling on droplet ac- 
tivation and condensational growth. We will refer to 
this as the “NO-QRAD’ ’ case, to distinguish it from 
the “base” case described in the previous section. Two 
direct effects are expected from omitting the effects of 
radiative heating on droplet condensational growth: de- 
creased condensational growth rates result in increased 
supersaturation and thereby increase droplet nuclea- 
tion, while increased values of 5 crit decrease droplet nu- 
cleation. The differences between the results are pro- 
nounced: compared to the base case, in the NO-QRAD 
case the total number of droplets at cloud top is more 
than doubled (Table 2). This increase is due to a 
greater concentration of droplets < 10 pm in radius. In 
the NO-QRAD case, the increase in droplet concen- 
tration results from a significant increase in the peak 
supersaturation. This result implies that, with regard to 
droplet nucleation, the effect of q TMi on condensational 
growth dominates the effect on S crjt . 

Other changes in cloud properties result from the 
increase in droplet concentration. Smaller average 
droplet size (as measured by r ert , the area- weighted 
mean radius, as well as r v ) leads to less efficient grow th 
of droplets by collection. Because droplet collection is 
an important sink for droplet number, less efficient col- 
lection magnifies the difference in droplet concentra- 
tion beyond that attributable to droplet activation alone. 
Decreased droplet collection also results in a nearly 
20% reduction in the peak drizzle rate. Lower drizzle 


rates allow the peak q x to increase by 24% in simulation 
NO_QRAD. The combination of smaller droplets and 
increased q x leads to a 38% enhancement in the cloud 
optical depth and relative increases of 12% in the 
broadband solar albedo and 45% in the peak longwave 
cooling rate. Greater cooling at cloud-top (an indirect 
result of omitting the effect of radiative transfer on 
droplet condensation and nucleation ) enhances vertical 
mixing, and the inversion height increases twice as 
much as it does in the base case. 

We conclude that in our model droplet radiative 
transfer has some very important effects on cloud prop- 
erties. But to the extent that droplets near the top of the 
cloud experience condensational growth for artificially 
long time periods (as discussed in section 3c), the 
model likely overestimates the effects of droplet radi- 
ative transfer. 

h. Infrared scattering 

Infrared scattering is typically ignored in cloud mod- 
els. To investigate its impact on our model results, we 
ran a simulation in which infrared scattering was set to 
zero ( NO_IR_SCAT). Without infrared scattering, in- 
frared photons have fewer opportunities to be absorbed 
by the cloud layer, and are more readily lost to the 
upper atmosphere. As seen in Table 2, ignoring infrared 
scattering leads to a 29% increase in the maximum 
longwave cooling. Greater cloud-top cooling increases 
the peak supersaturation sufficiently to activate smaller 
CN, thereby providing more droplets at cloud top. The 
increase in droplet concentrations leads to a small de- 
crease in r efT , which causes the drizzle rate to decrease 
slightly. Decreased drizzle and increased supersatura- 
tion allow the peak q x to increase by 10% over the base 
case. Optical depth and albedo are also increased. 

We conclude that infrared scattering by clouds is im- 
portant in our model and should require further inves- 
tigation with other models. 

c. Coalescence efficiencies 

The likelihood of droplets combining due to differ- 
ences in their fall speeds is represented by a matrix of 
gravitational collection efficiencies. The collection ef- 
ficiency between two droplets is the product of the like- 
lihood that they collide (the collision efficiency) with 
the likelihood that their collision results in their merg- 
ing into a single droplet (the coalescence efficiency). 
While coalescence efficiencies are often assigned a 
value of unity (e.g., Hall 1980; Flossman et al. 1985), 
in our model we employ the parameterization of coa- 
lescence efficiencies derived by Beard and Ochs 
( 1984). The smallest values of these efficiencies are 
for collisions between large droplets, where the coales- 
cence efficiency falls to a value of 0.5. To investigate 
the significance of these coalescence efficiencies, we 
ran a simulation (COAL = 1 ) in which all coalescence 
efficiencies were assigned a value of unity. 
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The differences between the COAL = 1 simulation 
and the base case are modest (Table 2). Near cloud 
base, where the drizzle flux is dominated by large 
drops produced from gravitational collection, the driz- 
zle flux for the COAL = 1 simulation is increased by 
— 10% over that in the base case. At the surface, the 
drizzle flux is increased by 65% over the base case. 
Although the peak in q ] is unchanged, the liquid water 
path (the vertical integral of q i) is reduced, resulting 
in a reduction of optical depth and small reduction of 
albedo. 

We conclude that realistic coalescence efficiencies 
should be used, although their effects on the micro- 
physics of stratocumulus clouds are relatively small in 
our model results. 

d. Collision efficiencies of small droplets 

There is some disagreement in the literature regard- 
ing the values of the collision efficiencies when both 
droplets are <30 pm in radius. Hall’s (1980) values, 
which are employed in our base case, were computed 
by Davis ( 1972) using Stokesian hydrodynamics. By 
approximating the effects of fluid inertia, Klett and Da- 
vis (1973) calculated increased collision efficiencies. 
For example, for a 20-pm radius droplet colliding with 
a 10-pm radius droplet, Klett and Davis (1973) ob- 
tained a collision efficiency that was twice the value 
calculated by Davis (1972) (13 and 6.8%, respec- 
tively). In simulation COLL1S<30 we investigated the 
sensitivity of our model results to these different col- 
lision efficiencies by using the values derived by Klett 
and Davis ( 1973 ). 

Simulation COLLIS<30 reduces the total concen- 
tration of droplets (Table 2), and there are more 
droplets > 20 pm in radius at all altitudes below 
cloud top than for the base case. The increase in 
larger droplets is reflected in increased drizzle fluxes. 
For example, at an altitude of 480 m, the drizzle flux 
increases 20% over that in the base case, while at the 
surface there is a 55% increase in the drizzle flux. 
Larger drizzle fluxes reduce q { . Lower q x and larger 
droplets (as measured by r cff ) combine to reduce the 
optical depth and the peak longwave cooling rate 
(both by 13%), and diminish the albedo (a relative 
decrease of 5% ). 

We conclude that a change in the collision efficien- 
cies of small droplets by a factor of two has noticeable 
effects on some of our model results. 

e. Condensation coefficient 

The condensation coefficient describes the fraction 
of water vapor molecules striking the surface of a drop- 
let that actually stick to the droplet. The value of the 
condensation coefficient is uncertain. For a quasi-qui- 
escent surface, an average value of 0.035 is suggested 
by Pruppacher and Klett ( 1978), while the values for 


a rapidly renewing water surface are larger by an order 
of magnitude and more. In our base simulation we use 
a value of 1.0 for the condensation coefficient. To in- 
vestigate the sensitivity of the model results to this pa- 
rameter, we ran a simulation in which a condensation 
coefficient of 0.035 was used (STICK=.035). 

As seen in Table 2, there is nearly a 50% increase in 
the peak value of the total droplet concentration in sim- 
ulation STICK=0.035 compared to the base case. Con- 
densational growth is inhibited in comparison to the 
base case, thereby limiting the rate at which small drop- 
lets grow large enough to be readily removed through 
collisions with larger drops. The drizzle flux is only 
slightly reduced (the greatest reduction, at cloud top, 
is <5%). Liquid water is increased at all altitudes 
above 400 m; there is a 10% increase in the maximum 
liquid water at cloud top. Decreased r cti and increased 
q x contribute to a 13% increase in both optical depth 
and the peak longwave cooling, and a 5% relative in- 
crease in albedo. 

We conclude that the value of the condensation co- 
efficient has a noticeable effect on cloud microstructure 
in our model. These effects are likely overestimated 
due to horizontal averaging, because droplets in the 
model are probably subjected to artificially long peri- 
ods of time in saturated conditions near cloud top (see 
section 3c ) . 

f Surface flux-profile relationships 

In our evaluation of surface fluxes we employ the 
dimensionless profiles of momentum and heat in the 
surface layer measured by Dyer ( 1974). An alternate 
set of profiles given by Businger et al. ( 1 97 1 ) was used 
by Duynkerke and Driedonks (1987). Besides sug- 
gesting slightly different dimensionless profiles, Bus- 
inger et al. ( 1971 ) also recommend a different value 
for von Karman’s constant (0.35, rather than 0.40, as 
usually assumed); furthermore, they concluded that the 
ratio of the eddy diffusivities for heat and momentum 
is >1. To investigate the sensitivity of our model re- 
sults to these differences, we ran a simulation 
(FLUX_PROFILES) that employs the dimensionless 
profiles of Businger et al., uses their recommended 
value of von Karman’s constant, and assumes a ratio 
of eddy diffusivity between heat (vapor also) and mo- 
mentum of 1.35 at all altitudes (<Th 1 = 1.35). 

Due to the smaller value of von Karman's constant, 
the surface shear production of E at 12 h is reduced by 
~30% in simulation FLUX_PROFILES, resulting in 
decreased E in the subcloud layer. But most of the dif- 
ferences between this simulation and our base case are 
due to the increased ratio between eddy diffusivities of 
momentum and vapor. The increased eddy diffusivities 
result in a greater flux convergence of q y in the upper 
region of the cloud layer, leading to increased super- 
saturation in this region (Table 2). The higher super- 
saturation nucleates more droplets and increases q x 
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above 600-m altitude, resulting in an enhancement of 
peak longwave cooling. In contrast to the other model 
simulations discussed above, where an increase in 
droplet concentration lowers drizzle fluxes and in- 
creases <?[, in this case an increase in q x leads to higher 
drizzle fluxes in the top 200 m of the cloud layer. Below 
600 m, where the flux divergence of q v is greater than 
in the base case, the supersaturation is lowered, reduc- 
ing q s and the concentration of droplets. These com- 
pensating effects result in essentially no change in liq- 
uid water path, optical depth, and albedo. 

We conclude that although the choice of surface flux 
profiles and eddy diffusivity ratios can affect mixing 
and cloud microstructure, cloud optical properties are 
not very sensitive to the values assumed for these pa- 
rameters. 

g. Probability distribution function for saturation 

Calculation of the buoyancy flux requires an evalu- 
ation of the fraction of air that is saturated at the bound- 
aries between grid layers (see appendix C). The frac- 
tional saturation is computed by integrating over all the 
saturated states in a semi -conservative thermodynamic 
phase space, which requires an assumed form of the 
distribution of states in the phase space. In our base 
case simulation, an exponential probability distribution 
function (PDF) is assumed. To explore the sensitivity 
of the model to the form of this function, we have run 
another simulation (PDF_HEAVISIDE), in which the 
PDF is represented by the Heaviside function. Instead 
of a continuous distribution of states, this step function 
describes a situation in which the air is either entirely 
saturated or entirely below saturation. This “all or 
nothing” condensation scheme is equivalent to that 
used by Duynkerke and Driedonks ( 1987 ) and Duyn- 
kerke ( 1988). 

The greatest impact of changing the PDF is greater 
decoupling between the cloud and subcloud layers. In- 
terestingly, N84 describes the turbulent structure of the 
atmospheric boundary layer during the period of ob- 
servations as decoupled. In comparison to our base 
case, in simulation PDF- HEAVISIDE q v is increased 
below 300-m altitude and slightly diminished above. 
Above 300 m, T is reduced by ~0. 1 K, and q { is slightly 
increased (Table 2). Because the buoyancy flux is in- 
creased in the lower regions of the cloud (where the 
fraction of condensed states in the base case is < I ), E 
is increased there in simulation PDF_HEAVISIDE. 
Below the supersaturated region of the main cloud 
layer (near 350-m altitude in this case), the fraction of 
condensed states in this simulation plummets to zero, 
resulting in a drastic minimum in turbulent fluxes that 
leads to a more pronounced minimum in E in this re- 
gion. Because vertical mixing between the main cloud 
and the subcloud layers is thereby limited, q y climbs 
enough in the subcloud layer to saturate one model 
layer in this simulation, at 300-m altitude. This feature 


could be interpreted to represent the cumulus clouds 
that N84 observed to rise into the thinning stratocu- 
mulus layer later in the afternoon. 

We conclude that the form of the probability distri- 
bution function used in the partial condensation scheme 
can have a noticeable effect on the vertical structure of 
the boundary layer. 

h . Treatment of horizontal divergence within the 

model domain 

Many models of the atmospheric boundary layer 
solve the dynamic equations in advective form and do 
not solve the continuity equation for air density (e.g., 
Bougeault 1985; Duynkerke and Driedonks 1987; 
Chen and Cotton 1987). When the vertical wind con- 
verges, air must either accumulate or “leak” out of the 
sides of the grid layers. When the ID dynamic equa- 
tions are solved in advective form, mass continuity is 
ignored, thereby implicitly leaking air from the sides 
to compensate for vertical mass convergence. Because 
the dynamic equations are solved in flux form in our 
model, in the base case air is explicitly leaked out of 
the sides, thereby keeping a constant air mass column 
in the model domain [ the second term in Eqs. ( 1 ), ( 2 ), 

( 10), ( 1 1 ), and ( 13)]. In a further simulation (PRES- 
SURE-RISE), we investigated the effect of allowing 
the air mass in a column to accumulate in the model 
domain by omitting this second term. For this sensitiv- 
ity test, pressure was recalculated at all levels after each 
time step. The pressure at the top of the model was 
fixed at its initial value, which implies that horizontal 
convergence above the model domain maintained the 
mass of the air column. Instead of balancing the fluxes 
of particles, vapor, and potential temperature out of the 
bottom of the highest model layer, advective fluxes into 
the top of the model were calculated by assuming con- 
stant vertical gradients of concentrations. 

In 12 h of the PRESSURE-RISE simulation, the 
surface pressure increased by 14 mb. Such a rate of 
pressure increase cannot be sustained for much longer 
than a day before the surface pressures become un- 
realistic (the leaking and accumulating treatments 
could be combined to offset this problem over longer 
simulations ). The changes in pressure and density raise 
the top of the model by 107 m from its initial altitude 
of I km. The increase in pressure does not affect 0 (at 
constant sigma levels) but it does increase T. Because 
q x is not affected by the increase in pressure, the rela- 
tive humidity decreases. The reduced relative humidity 
decreases cloud depth (Table 2) and evaporates the 
drizzle before it reaches the surface. The liquid water 
path is reduced significantly, leading to reductions in 
optical depth and the peak longwave cooling rate. The 
resulting decrease in albedo is small. 

We conclude that the method of treatment of hori- 
zontal divergence modestly affects our model results, 
and it has a pronounced effect on drizzle reaching the 
surface. 
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Table* 5. Variations in input CN si/e distributions. The volume-weighted mean radius (r v ) is determined by and a. 


Model simulation 

r„ 

(pm) 

( T 

t\ 

(pm) 

Initial 

concentration 
(cm ■') 

Production 

rate 

(cm * s 1 ) 

BASE-CASE 

0.05 

2.5 

0.62 

1000 

0.018 

CN_SIG2_R.()5 

0.05 

2.0 

0.21 

1600 

0.028 

CN_SIG2_R. 15 

0.15 

2.0 

0.65 

500 

0.005 

CN_S1G1.2_R.12 

0.12 

1.2 

0.15 

500 

0.005 

CN_S1G1.2_R.56 

0.56 

1.2 

0.62 

120 

0.0055 


5. Sensitivity of model to environmental parameters 

The environmental parameters measured by N84 are 
nearly sufficient to specify all of the boundary condi- 
tions required to run our base case model simulation 
(section 3a). However, the values of a few parameters 
had to be assumed. Notable among these is the char- 
acterization of the CCN in the atmosphere during the 
period of observations. Also, the divergence rate of the 
horizontal wind velocity, the sea surface temperature, 
and the history of the air mass were not reported by 
N84. There is also uncertainty in the downwelling flux 
of longwave radiation. In this section we investigate 
the sensitivity of the model results to the assumed val- 
ues of these environmental parameters. 

a. Characterization of CCN 

Nicholls ( 1984 ) did not measure droplets < 1 // m in 
radius, therefore nothing is known about the CCN upon 
which the observed cloud droplets formed. As de- 
scribed above, for the input CN distribution in our base 
case simulation we matched the shape of CCN activa- 
tion spectra measured underneath stratus clouds off the 
California coast and adjusted the total number and pro- 
duction rate to approximately match the average total 
droplet concentration measured by N84. To examine 
the effects of variations in the input CN distribution on 
the model results, we ran five further simulations: one 
to evaluate the role of sea salt particles, two to address 
the sensitivity of the results to small variations in the 
slope of the CCN spectrum, and two to investigate the 
effects of nearly monodisperse input CN distributions. 
Some caution is warranted in interpreting the signifi- 
cance of the results of these sensitivity tests, given the 
importance of supersaturations on droplet nucleation 
and the likely shortcomings of the horizontally aver- 
aged supersaturation field predicted by the model. 

Because the spectral width of the CN distribution 
used in our base case simulation is so broad (a ~ 2.5 ), 
it represents not only the small sulfate particles (r < 1 
//m) that dominate the number concentration of parti- 
cles in the marine atmosphere but also the “giant' 1 nu- 
clei (r > 1 jum) consisting of sea salt particles. This 
was verified by comparing our CN distributions with 
sea salt measurements by Woodcock ( 1953 ) taken at 
wind speeds corresponding to those observed by N84. 


Due to interest in possible effects of giant nuclei on 
precipitation (e.g.. Beard and Ochs 1993), we evalu- 
ated their role under the conditions observed by N84 
by running a model simulation ( NO_SEA_SALT ) in 
which there were no input CN with radius > 1 /zm. In 
our base case simulation there are initially 1 cm * of 
CN with radii > 1 pm (out of 1000 total), and giant 
nuclei compose 35% of the total mass. 

The effects of eliminating the source of giant nuclei 
are modest (Table 2). The number of droplets < 5 /xm 
in radius increases, leading to an increase in the total 
concentration of droplets ( /V drops ) . Although the peak 
supersaturation at cloud top is the same as in the base 
case, supersaturation below cloud top is slightly in- 
creased in the NO_SEA_SALT simulation. This in- 
crease in supersaturation leads to a slight increase in 
r ett that gives rise to an increased drizzle flux. The re- 
sulting relationship is a counter example to the con- 
ventional wisdom that increased /V arnps leads to de- 
creased drizzle fluxes (Albrecht 1989; Baker and 
Chari son 1990). The difference between our results in 
this case and conventional wisdom is because in our 
model the increase in droplet concentration is limited 
to droplets that are so small that they have little effect 
on the cloud microphysics. We conclude that under 
conditions of relatively high droplet concentrations, gi- 
ant nuclei play a minor role in the microphysics of our 
model. 

The input CN distributions for the four remaining 
CN sensitivity tests are given in Table 3. In simulations 
CN_SIG2_R.05 and CN_SIG2_R.15, the geometric 
standard deviation of the CN distribution (a) is re- 
duced from 2.5 to 2; in simulation CN_SIG2_R.05, r n 
is unchanged from the base case, while in simulation 
CN_SIG2_R. 15, the geometric mean volume radius 
( r v ) is unchanged. The initial number and production 
rate of CN were chosen to produce a below-cloud CCN 
activation spectra at 12 h that intersects that of the base 
case at its peak supersaturation of —0.04% ( seen in Fig. 
11), since such a constraint could be expected to result 
in similar /V drops . The issue of interest in these two sim- 
ulations is to see whether or not small changes in the 
slope of the CCN activation spectrum affect the model 
results. 

The activation spectrum (for supersaturation be- 
tween 0.01 and 0.1%) in simulation CN_SIG2_R.05 
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is steeper than that for the base case (Fig. II). The 
increased slope results in increased supersaturation at 
all altitudes within the cloud layer, which produces 
higher above 600-m altitude. However, the higher 
supersaturations also produce an increase in r L , ff . Bigger 
droplets produce greater drizzle fluxes than in the base 
case, as well as reductions in optical depth, peak long- 
wave cooling, and albedo. This simulation is another 
counter example to the conventional wisdom that in- 
creased /Vjrops leads to reduced drizzle fluxes. 

In contrast to the previous simulation, in simulation 
CN_SIG2_R.15 the activation spectrum for supersat- 
urations between 0.01 and 0.1% is less steep than for 
the base case (Fig. 11). The decreased slope results in 
lower supersaturations, and A drops is increased at all al- 
titudes. The increase in A drups results in slightly reduced 
r eff in the upper region of the cloud layer. Although the 
drizzle fluxes are changed only slightly, the cloud layer 
is deeper than in the base case, giving rise to increased 
cloud water. The increased q } leads to a greater optical 
depth and albedo. We conclude that the small variations 
in the distribution of CN have a modest effect on our 
model results. 

The input size distributions for simulations 
CN_SIG1.2_R.12 and CN_SIG1.2_R.56 are nearly 
monodisperse (a = 1.2). The particle size distribution 
in CN_SIG1.2_R.12 is similar to surface measure- 
ments of sulfate over the remote, cloud-free Pacific 
Ocean (Clarke et al. 1987). The size distribution of 
CN in CN_S1G1.2_R.56 represents oversimplification 
of the nucleation process in that either all or none of 
the CN are activated. The peak supersaturation from 
the base case ( —0.04%) corresponds to the maximum 
slope in the activation spectrum of CN_SIG 1 .2_R. 12 
(Fig. II), while all of the CN are activated at that su- 
persaturation for the distribution in CN_SIG1 .2_R.56. 

In simulation CN_SIG 1 .2_R. 1 2, the supersaturation 
is increased at all levels in the cloud layer, but nucle- 
ation is limited to cloud top, and is increased only 
in the upper half of the cloud layer. At cloud top the 
dramatic increase in A drops is caused by a new peak in 
the size distribution for droplet radius between 2 and 3 
^m. But these droplets are so small that their impact 
on cloud microphysics is overwhelmed by the effect on 
condensational growth of increased supersaturations 
throughout the cloud layer. This enhanced condensa- 
tional growth produces increased r e(t and drizzle fluxes. 
Increased drizzle decreases q \ , this, in combination 
with increased r cff , results in a 15% decrease in optical 
depth and a relative decrease in albedo of 13%. This 
simulation serves as yet another counter example to the 
expected effects of increased A drops on drizzle, and 
albedo. 

In simulation CN_SIG1.2_R.56, the peak of super- 
saturation is drastically reduced (by a factor of three) 
from that in the base case, but, because the CN are 
activated at such a low supersaturation ( below 0.01%), 
A^s is enhanced over that in the base case. In this 
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Fig. 1 1. Comparison of modeled cumulative CCN activation spec- 
tra at 12 h and 150-m altitude for base case (solid). CN_S1G2_R.()5 
(dotted). CN-SIG2—R. 15 (dashed), CN-SIGl .2^R.12 (dot dashed), 
and CN-S1GI.2-R 56 is tri-dot -dashed. 


simulation the increase in A drops results in a decreased 
r efT , which leads to decreased drizzle fluxes, resulting 
in increased peak q } . Smaller droplets and increased q\ 
reinforce an increase in optical depth and albedo. The 
increase in optical depth leads to an increased peak 
longwave cooling, which results in a deeper cloud that 
entrains more inversion air. We conclude that treating 
the marine aerosol as nearly monodisperse has a sig- 
nificant impact on the model results. 

b. Downwelling longwave radiation 

As is typical of marine stratocumulus, in the obser- 
vations reported by N84, as well as in our model sim- 
ulations, the buoyancy flux provided by cloud-top ra- 
diative cooling drives most of the vertical mixing in the 
boundary layer. Cloud-top radiative cooling depends 
on the downward flux of longwave radiation, which is 
specified in our base case simulation (at 270 W m" 2 ) 
to match the radiative calculations presented in N84. 
N84 also measured the downwelling longwave radia- 
tive flux, which had an average value 15 W m 2 greater 
than their calculations. In their study of five further 
cases of stratocumulus over the North Sea, Nicholls and 
Leighton ( 1986 ) found a similar error of comparable 
magnitude in all of their comparisons of radiative cal- 
culations and measurements for downwelling radiative 
flux above cloud top, which they ascribed to a system- 
atic measurement error. Because this flux is so impor- 
tant to the cloud-top radiative cooling rate, we ran a 
model simulation (LW-FLUX) in which the down- 
welling flux at 9(X)-m altitude was the same as that 
measured in N84 ( at 285 W m 2 ) . This increase in the 
downwelling flux above cloud top was produced by 
increasing the temperature of the blackbody radiative 
source at the top of layer 0. 

After 12 h of this simulation, the peak longwave 
cooling rate is 22% less than in the base case, which 
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reduces vertical mixing. The in-cloud peak of E is re- 
duced by 20% from the base case, and the inversion 
height increases by only 10 m in 12 h. Less mixing and 
less cooling produce a 14% decrease in the peak su- 
persaturation, though the peak in r/, is only slightly re- 
duced. Because of reduced mixing, the drizzle fluxes 
are diminished at all altitudes. The changes in mixing 
and drizzle bring the model results in closer agreement 
to the observations. We conclude that a 15 W m 2 in- 
crease in the downwelling longwave radiative flux has 
an appreciable effect on cloud dynamics, but the effect 
on optical properties is modest. 

Chen and Cotton ( 1987) performed a similar sensi- 
tivity test of their model results to the downwelling 
longwave radiative flux by prescribing an upper-level 
cloud deck above the boundary layer ( their experiment 
HIC1 ). The increase in downwelling longwave radia- 
tion in their sensitivity test was greater than in ours ( it 
reduced their peak longwave cooling rate by >50%), 
and resulted in a greater decrease in the peak value of 
q\ (it fell by nearly 30% in their test). 

c. Divergence rate of horizontal wind velocity 

Because the large-scale divergence rate of horizontal 
wind velocity (div) represents a small difference be- 
tween large terms, it is difficult to measure reliably. In 
the observations of N84 it was not reported. In our base 
case simulation we use a value of 2.5 X 10 6 s 1 . By 
causing subsidence of (warm and dry) inversion air, it 
affects the dynamics of boundary layer mixing by off- 
setting cloud-top longwave cooling. To investigate the 
effects of variations in the divergence rate, we have run 
two further simulations: DIV0, in which no divergence 
is imposed, and DIVX2, in which div is doubled from 
the value used in the base case. 

In the base case simulation, subsidence produces 3. 1 
K h 1 of heating at cloud top at 12 h, which is more 
than one-third the longwave cooling rate (Table 2). In 
DIV0 this subsidence heating no longer operates, 
which allows the inversion height to increase by twice 
that in the base case and results in a cloud layer over 
100 m deeper. Although the cloud layer is deeper, in- 
creased drizzle at all altitudes prevents any increase in 
the peak value of q K . But the increased cloud depth does 
increase the liquid water path, resulting in a 11% in- 
crease in optical depth and a relative increase in albedo 
of 4%. 

Doubling the divergence rate to 5 X 10 -6 s” 1 in sim- 
ulation DIVx2 leads to a more sharply defined inver- 
sion at 12 h, in which the increased subsidence pro- 
duces 7.2 K h 1 of heating at cloud top. This heating 
nearly offsets the cloud-top longwave cooling, al- 
though the latter extends down a few more layers, 
where it is not offset by subsidence. Less net cloud-top 
cooling results in a cloud layer that is 80 m shallower 
than in the base case. Instead of the inversion height 
increasing, it decreases by 20 m below its initial value 


in simulation DIVX2. With less net cloud-top cooling, 
the peak supersaturation is reduced; this, in combina- 
tion with the shallower cloud layer, reduces the peak t/i 
by 12%. The total concentration of droplets increases 
because drizzle rates are significantly reduced at all al- 
titudes. The net of these changes is to produce a 15% 
reduction in optical depth and a relative albedo reduc- 
tion of 5%. We conclude that variations in the diver- 
gence rate can cause significant variations in cloud 
properties. 

Chen and Cotton ( 1987 ) found a similar sensitivity 
of their model results to the divergence rate of hori- 
zontal wind velocity. In their simulations of a shal- 
lower, shear-driven stratocumulus-topped marine bound- 
ary layer, they found that the peak value of q x decreased 
in response to an increased divergence rate. 

d. Sea surface temperature and airmass history 

Because sea surface temperatures (SST) were not 
reported by N84, for the base case simulation we used 
the climatological average July value for the measure- 
ment area. To investigate the sensitivity of our model 
results to a small difference in SST, we ran a further 
simulation in which SST was increased by 1 K. An- 
other unknown in comparing our model results with the 
N84 observations is the history of the air mass. Al- 
though the measurements presented in N84 character- 
ized the state of the boundary layer over a period of 
~4 h, they were still only a snapshot in time. To 
crudely represent an alternative history of the air mass 
from that represented by our base case simulation, we 
ran a further simulation in which the initial tempera- 
tures throughout the boundary layer were increased by 
I K. This alternative history is logistically designed to 
bring the predicted boundary layer temperatures closer 
to the observations; it is physically justifiable because 
the air mass may have been influenced by passing over 
Scotland before it reached the North Sea. 

In simulation SST+1 the SST is increased by 1 K. 
After 12 h, SST+1 results in a surface buoyancy flux 
of 19 W m" 2 , compared to 14 W m 2 in the base case. 
The enhanced surface buoyancy produces turbulent 
fluxes near the surface that are in the upper range of 
the observed fluxes. In contrast, the fluxes near the sur- 
face in the base case split the differences between the 
observed values (Fig. 8). The increased surface source 
of heat (both sensible and latent) produces a boundary 
layer at 12 h that was ~0.5 K wanner than in the base 
case; this matches the temperature and vapor measure- 
ments more closely, but is still colder ( by —0.5 K ) than 
the measurements in the cloud layer. In contrast, 
SST+1 produces an upwelling longwave radiative flux 
near the surface that is further from the radiative cal- 
culations of N84 than our base case simulation, in 
which the values are already slightly higher than N84 
calculated (Fig. 6). Other changes from the base case 
include a 10% increase in the peak value of the tur- 
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bulent kinetic energy in cloud, an inversion height in- 
crease 20 m greater than in the base case, and slightly 
increased drizzle fluxes. Although Table 2 indicates 
that the peak longwave cooling rate is 22% less than 
the base case, this is a transient effect due to the en- 
trainment of a new layer just before the output at 12 h; 
it does not represent a sustained difference from the 
base case. 

In testing the sensitivity of their model results to an 
increased SST ( which they increased by 7 K, in con- 
trast to our increase of IK), Chen and Cotton ( 1987) 
also found that entrainment of inversion air was en- 
hanced. However, in contrast to our results, they found 
that the peak value of q\ was reduced, and the cloud 
base descended, when the SST was increased ( in our 
results neither changed in response to the smaller in- 
crease in SST). 

In simulation BLT+ 1, the initial temperatures in the 
boundary layer are increased by 1 K, but the SST is the 
same as the base case. This reduces the surface bound- 
ary flux to 9 W m 2 at 12 h (compared to 14 W m 2 
in the base case). The boundary layer temperatures at 
12 h in simulation BLT-Fl are —0.1 K higher than in 
simulation SST+I. Although there is a slight increase 
in the droplet concentration, as in simulation SST+1, 
the changes in cloud optical properties are negligible. 
We conclude that small variations in SST and initial air 
temperature in the boundary layer produce small 
changes in the model results. 

6. Summary 

In this paper we have described a I D model for the 
stratocumulus-topped marine boundary layer, which 
includes aerosol and cloud microphysics, radiative 
transfer, and vertical mixing. We have compared the 
outputs of the model with measurements and have 
tested the sensitivity of the model to physical assump- 
tions and environmental uncertainties. 

Although the model includes a detailed treatment of 
aerosol and cloud microphysics and radiative transfer, 
the representation of air motions is highly simplified. 
The most significant shortcoming of 1 D models such 
as this is that grid layers represent averages over up- 
drafts and downdrafts. In the upper region of the cloud 
this model simplification probably results in droplets 
being exposed to supersaturated conditions for artifi- 
cially long time periods. In the lower region of the mod- 
eled cloud, on the other hand, there is probably exces- 
sive evaporation of cloud droplets, which likely results 
in overprediction of mixing between the cloud and the 
subcloud layers. Another likely consequence of aver- 
aging over updrafts and downdrafts is that peak super- 
saturations are underpredicted by the model. Therefore, 
the model likely underpredicts the number of CN that 
serve as CCN. 

With these caveats in mind, the following results 
have been obtained from the model simulations de- 
scribed in this paper. 


• The model predictions of thermodynamics, micro- 
physical properties, and radiative fluxes are generally 
in reasonable agreement with measurements in an 
~500-m thick, summertime marine stratocumulus 
cloud layer. However, the model predictions of turbu- 
lent fluxes between the cloud and subcloud layers ex- 
ceed the measurements. 

• The variations in cloud properties between thick 
stratus measured under gale conditions and high, thin 
stratocumulus layer are generally reproduced by the 
model, although it underpredicts the entrainment of 
overlying air at cloud top under gale conditions. Al- 
though the drizzle fluxes predicted by the model are 
within the measurement uncertainties, the maxima are 
consistently overpredicted. 

• Most observations of stratiform clouds show that 
the total concentration of cloud drops (N drops ) is not 
strongly dependent on height. However, the modeled 
profile of 7V dmps is not necessarily constant with height 
in the cloud layer. Changes in /V drups with height are 
sensitive to the lower droplet size cutoff and the pop- 
ulation of large unactivated haze particles. For a lower 
droplet size cutoff of 2.2 pm, the profile of droplet 
concentration correlates with the supersaturation pro- 
file in the cloud layer. Particles with radius ~1 pm 
can contribute to /V dmps , but because they are so small 
they have little effect on the microphysics of the cloud 
layer. 

• While it is often the case that increases in N dri>ps 
lead to decreased drizzle fluxes, and therefore increases 
in cloud water, this is not always the case in our model 
simulations. This is because in our model increased 
concentrations of small droplets (r < 5 pm) can en- 
hance /V dn , ps while having little effect on cloud micro- 
physical processes. When this happens, the radius of 
the mean droplet volume (r v ) is no longer a useful 
indicator of average droplet size (i.e., increased r v is 
not associated with increased drizzle fluxes), because 
it is so sensitive to the concentrations of small droplets. 
The area- weighted droplet radius ( r cii ) is a more mean- 
ingful indicator of the average size of cloud droplets 
because, in our model simulations, it was always pos- 
itively correlated with the drizzle flux. 

• The effect of radiative transfer on droplet conden- 
sational growth leads to significant changes in the 
structure of the modeled cloud layer. Droplet radiative 
cooling leads to a lower supersaturation and less drop- 
let nucleation near cloud top. This effect is generally 
ignored in cloud models; our model results suggest that 
it may be important when cloud-top radiative cooling 
is significant. We recommend that this effect be ex- 
plored further with other models. 

• Infrared scattering, which is also generally ignored 
in cloud models, has noticeable effects on the structure 
of the modeled cloud layer. Cloud-top radiative cooling 
decreases when infrared scattering is included in the 
model. 
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• A factor of 2 change in the collision efficiencies 
between droplets with radii <30 pm has a noticeable 
effect on the properties of the modeled cloud layer. 

• The value of the condensation coefficient for drop- 
let growth has a noticeable effect on the properties of 
the modeled cloud layer; a smaller value produces 
smaller droplets and increased cloud water. 

• Small differences in the characterization of the in- 
put particle distributions influence the model results 
only slightly. However, the assumption of a nearly 
monodisperse CN distribution dramatically alters the 
properties of the modeled cloud layer. 

• The divergence rate of the horizontal wind velocity 
has a noticeable effect on the modeled cloud structure. 
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APPENDIX A 

Symbols 

A k , A s Kelvin and solute factors in the conden- 
sational growth kernel 
A v Avogadro’s number 

B mean radiative intensity 

c d , c h surface drag and heat transfer coefficients 

c p specific heat at constant pressure of dry air 

C concentration of droplets in a size bin 

c'i, , c 2t constants in turbulence model [values 
given in Duynkerke ( 1988 ) | 

D water vapor diffusion coefficient 

div divergence rate of mean horizontal wind 

velocity 

E turbulent kinetic energy 

/ Coriolis parameter 

F v , F x vapor and thermal ventilation factors 

g gravitational acceleration 

G water vapor concentration 

g r droplet condensational growth rate 

J Planck function 

K m gradient diffusion coefficient for momen- 

tum 

k von Karman’s constant 

k dbs radiative absorption coefficient for a par- 

ticle size bin 

K b , K c thermal (Brownian) and total coagulation 
kernels 

K y thermal conductivity of air 

L Monin-Obukhov length 

L v latent heat of water vaporization 


M 

mean wind speed 

m s , 

masses of solute and water, respectively, 
in a droplet 

M s , M w 

molecular weights of solute and water 

n v 

ambient vapor pressure 

^drops 

total droplet concentration 

rtvup 

saturation vapor pressure of liquid water 

p 

production rate of E 

</rad 

droplet radiative heating rate 


saturation vapor mixing ratio 

<7t 

total water mixing ratio 

</v 

liquid and water vapor mixing ratios 

r, r' 

particle radii 

R 

universal gas constant 

Cut 

critical radius for droplet activation 

h 

radius of mean droplet volume 

Cn 

effective (area-weighted) droplet radius 

Cm < 

geometric mean number and mean vol- 
ume radii, respectively, of lognormal 
CN distribution 

Rn 

particle loss rate 

s 

water vapor supersaturation 

*^crit.Ki ^ri! 

critical supersaturations for droplet acti- 
vation, ignoring and including radiative 
effect 

s n 

source of particles 

T, T 0 

time-dependent and initial air tempera- 
tures 

t 

time 

M, U 

horizontal wind components 


surface wind stress 

'V U g 

components of geostrophic wind 

Clep.p ^dcp.v 

deposition velocities for particles and wa- 
ter vapor 

t^dcp.W 

surface transfer velocity for sensible heat 

t’f 

particle fall speed 

W 

vertical wind speed 

w* 

convective velocity scale 

z 

altitude above ocean surface 

Z{) 

surface aerodynamic roughness height 

Zi 

altitude at which surface values of E and 
e are evaluated ( 1.5 m) 

Zh 

altitude at which stress falls to 5% of its 
surface value 

c 

dissipation rate of E 

<$> 111 , <t> h 

dimensionless gradients of momentum 
and potential temperature 


practical osmotic coefficient of dissolved 
CCN 

V, \ 

radiative frequency and wavelength 


dissociativity of dissolved CCN 

0 

potential temperature 

0* 

equivalent potential temperature [= 0 + 
(L v 0/c p T)q v ] 

0 { 

liquid potential temperature [ = 0 - ( L v 0 / 
c p T)q x ] 

^v() 

time-dependent and initial virtual poten- 
tial temperatures 
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P . Pw 

densities of air and liquid water 

a 

geometric standard deviation of lognor- 
mal CN distribution 

•> 

turbulent Prandtl numbers for particles 
and heat 

a B , & f 

turbulent Prandtl numbers for E and t 
[ values given in Duynkerke ( 1 988 ) ] 

(7 w 

surface tension of water 

<A,n. 'Ah 

integrated dimensionless profiles of mo- 
mentum and temperature 

APPENDIX B 

Buoyancy Flux 


The weighting R s is determined through a partial con- 
densation scheme, in which an estimate is made of the 
fraction of saturated states in the phase space of the 
semiconservative variables 0\ and q { . The formulation 
assumes that all liquid water is instantly available to 
maintain saturation, which is not necessarily the case 
in unsaturated regions. The method has been described 
by Bougeault ( 1981 ). Here R s is found by integrating 
a normalized dimensionless variable in — q y phase 
space that describes the distribution of states; its stan- 
dard deviation is given by 

+ aW? -2a^[0\) u \ (B7a) 


The buoyancy flux is an important term in the E 
equation for the stratocumulus-topped marine bound- 
ary layer. Because 0 V is not conserved under saturated 
conditions, gradient transfer is only appropriate under 
unsaturated conditions. Here we describe how the sat- 
urated and unsaturated fluxes are combined. 

Applying a Reynolds average to the definition of 0 V 
yields its turbulent flux: 


where 


a = 1 + 


Uq si \ 
c p R v Tl) 


and 


_ T Lygsl 
_ (9 RJ} ' 


( B7b) 


To evaluate Eq. ( B7 ) we use the relations from the 2.5 
order closure model of Mellor and Yamada (1982), 
which simplify to 


w'9' y = (1 + 0.6k/ v - q ] )w , O t 

T 0.61#vc f q[ - Ow f ql • (Bl) 


, a 2 B 2 lK(dq t dO t \ 2 
4,l2E°»\dz a, dz) 


( B8) 


Under saturated conditions, the Clausius-Clapeyron 
relation and the Boussinesq approximation lead to 


w'q [) sat = a w'0\ where a 


T d % _ L y q s 
0 d T ~ RJO ' 

(B2) 


and q s is the saturation specific humidity. Applying a 
Reynolds average to a linearized form of <9 e : 

OE 

w ' 9 ' = vv ’0' c w ’q[. ( B3 ) 

Tc P 

Equations (B2 ) and (B3) yield 
w T ql) M = ^w r 5l, where j8 = I +^a. (B4) 

fj I C p 

Substituting Eq. (B4) back into Eq. (B3) yields 

ifTL = WWjp. ( B5 ) 

Finally, the turbulent flux of q x is just the difference 
between the turbulent fluxes of total water (q t ) and <? v . 
The saturated buoyancy flux is evaluated through Eqs. 
(B3)-(B5), in which the turbulent fluxes of the semi- 
conservative variables 0 C and q y are evaluated through 
gradient transfer, while the unsaturated buoyancy flux 
uses gradient transfer directly. 

The unsaturated and saturated buoyancy fluxes are 
combined through a weighting factor, R s , which rep- 
resents the fraction of saturated air at a given level: 

(B6) 


where / = c'J* E* /2 fe is the master mixing length, and 
B 2 — (2/3) (4/c^) V4 . For the probability distribution 
function, we use the exponential form given by Bou- 
geault ( 1980). As discussed by Bougeault ( 1981 ), the 
exponential form differs from the Gaussian form only 
at extreme values (nearly 0 or 1 ) of the normalized 
dimensionless variable in 9 { - q { phase space, and then 
only in a relative sense, as the absolute value of R s is 
very nearly 0 or 1 either way. Such subtle differences 
in the probability distribution function make little dif- 
ference to our model results. 


APPENDIX C 


Surface Similarity Boundary Conditions 


Surface fluxes depend upon the drag coefficient c d 
= u\l M 2 , in whic h the surface stress is defined by 
u% = {{w'u'jl T {W'v')l] u2 . The mean wind speed 
in the lowest model layer is calculated by integrating a 
dimensionless stability-dependent wind speed profile: 


<f>m(zfL) = 


kz dM 


(Cl) 


where k is von Karman’s constant and L the Monin- 
Obukhov length: 


u * 9 Vii 

kgWOl) s’ 


(C2) 


which is the ratio of surface shear production of E to 
the surface buoyancy flux. The integration of the wind 


w’0’ v = R s w'0;.) sat + ( 1 - R s )W* )unsat • 
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profile from z 0 (the aerodynamic roughness length) to 
z yields 



thereby allowing the drag coefficient to be evaluated 
from c‘ d ~ ip m 2 . Similarly, a transfer coefficient for 
transport of heat and vapor across the turbulent (con- 
stant flux ) sublayer is evaluated from c h - t/'m 1 iK 1 ■> in 
which is an integrated dimensionless temperature 
profile. For the dimensionless profiles of wind speed 
and temperature, we follow Duynkerke ( 1988 ) and use 
the measurements of Dyer ( 1974). To evaluate the in- 
tegrated profiles, we use the numerical method devised 
by Benoit ( 1977 ) for the unstable surface layer, which 
avoids numerical difficulties when approaching the 
neutral limit. 

The aerodynamic roughness length is evaluated at 
the beginning of each time step from Charnock's 
( 1955 ) relation: 



in which b c 1 = 0.015. To iteratively calculate the sur- 
f ace fl uxes within each time step, we first evaluate 
(w'Ol ) s using the previous value of c\, and the current 
values of A/, #(z s ), and q K ( z s ) ■ Then L is evaluated, 
which allows a reevaluation of c\ t and c h , followed by 
a rec alculation of deposition velocities. Finally 
(w , '#' ) s is reevaluated, and compared to its most re- 
cently calculated value . This loop iterates until the rel- 
ative change in (w f 0[. \ falls below 10 6 . 
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